diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index bc25b97012..9890f63dfd 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -462,6 +462,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])())); diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C index 5412ef5b1c..ac7eac71c1 100644 --- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C +++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C @@ -439,6 +439,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C index 207d3d9cbc..68ea711a06 100644 --- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C +++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C @@ -332,6 +332,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" bool overwrite = args.options().found("overwrite"); diff --git a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C index 33d9ba71af..0fcea92bfe 100644 --- a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C +++ b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C @@ -56,6 +56,7 @@ int main(int argc, char *argv[]) argList::validArgs.append("cellSet"); # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" pointMesh pMesh(mesh); diff --git a/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C b/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C index f4336bce92..6f476206fb 100644 --- a/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C +++ b/applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C @@ -54,6 +54,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" word patchName(args.additionalArgs()[0]); diff --git a/applications/utilities/mesh/advanced/removeFaces/removeFaces.C b/applications/utilities/mesh/advanced/removeFaces/removeFaces.C index 5cf8473b5d..b1b5695c70 100644 --- a/applications/utilities/mesh/advanced/removeFaces/removeFaces.C +++ b/applications/utilities/mesh/advanced/removeFaces/removeFaces.C @@ -53,6 +53,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" bool overwrite = args.options().found("overwrite"); diff --git a/applications/utilities/mesh/advanced/splitCells/splitCells.C b/applications/utilities/mesh/advanced/splitCells/splitCells.C index b5a67074e6..be28bf6893 100644 --- a/applications/utilities/mesh/advanced/splitCells/splitCells.C +++ b/applications/utilities/mesh/advanced/splitCells/splitCells.C @@ -532,6 +532,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); diff --git a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C index e982206e61..ffa8314952 100644 --- a/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C +++ b/applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C @@ -47,6 +47,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); diff --git a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C index 1d2eaae78b..59423a9d04 100644 --- a/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C +++ b/applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C @@ -346,6 +346,7 @@ int main(int argc, char *argv[]) # include "addTimeOptions.H" # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); bool patchFaces = args.options().found("patchFaces"); bool doCell = args.options().found("cell"); diff --git a/applications/utilities/mesh/generation/extrude2DMesh/doExtrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/doExtrude2DMesh.C index d974dc8da2..30a8bce0fc 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/doExtrude2DMesh.C +++ b/applications/utilities/mesh/generation/extrude2DMesh/doExtrude2DMesh.C @@ -61,6 +61,7 @@ int main(int argc, char *argv[]) argList::validOptions.insert("overwrite", ""); # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" scalar thickness(readScalar(IStringStream(args.additionalArgs()[0])())); diff --git a/applications/utilities/mesh/manipulation/attachMesh/attachMesh.C b/applications/utilities/mesh/manipulation/attachMesh/attachMesh.C index b808dd7521..2b0f894e98 100644 --- a/applications/utilities/mesh/manipulation/attachMesh/attachMesh.C +++ b/applications/utilities/mesh/manipulation/attachMesh/attachMesh.C @@ -46,6 +46,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" bool overwrite = args.options().found("overwrite"); diff --git a/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C b/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C index 3ea6a2a02b..f7cea3d4a0 100644 --- a/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C +++ b/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C @@ -75,6 +75,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" Info<< "Mesh read in = " diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C index 2f1fd57634..882b7461e5 100644 --- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C +++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C @@ -58,6 +58,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" const polyBoundaryMesh& patches = mesh.boundaryMesh(); diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C index 9ba825083f..a449661982 100644 --- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C +++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C @@ -308,6 +308,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); const bool overwrite = args.options().found("overwrite"); diff --git a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C index 8ec01e2920..cbd2d9e419 100644 --- a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C +++ b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C @@ -160,6 +160,7 @@ int main(int argc, char *argv[]) argList::validOptions.insert("overwrite", ""); # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" bool split = args.options().found("split"); diff --git a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C index 33bb0e00d9..8cde78ceb4 100644 --- a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C +++ b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C @@ -298,6 +298,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" printEdgeStats(mesh); diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index be3a2debb7..c7e85f0a6e 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -374,6 +374,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); // Get times list instantList Times = runTime.times(); diff --git a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C index d902303bca..38a09fe2e1 100644 --- a/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C +++ b/applications/utilities/mesh/manipulation/splitMesh/splitMesh.C @@ -122,6 +122,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createPolyMesh.H" word setName(args.additionalArgs()[0]); diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C index 511141ac66..eb116fc231 100644 --- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C +++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C @@ -1122,6 +1122,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" word blockedFacesName; diff --git a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C index e3c558ce50..4452bf46ca 100644 --- a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C +++ b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C @@ -135,6 +135,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C index d49fd4c6b5..b245f0bd0b 100644 --- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C +++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C @@ -157,6 +157,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" + runTime.functionObjects().off(); # include "createMesh.H" word setName(args.additionalArgs()[0]); diff --git a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C index 50bf81035f..7bff5f4c23 100644 --- a/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C +++ b/applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C @@ -100,13 +100,23 @@ int main(int argc, char *argv[]) if (dict.found(entryNames[0])) { - const entry* entPtr = &dict.lookupEntry(entryNames[0]); + const entry* entPtr = &dict.lookupEntry + ( + entryNames[0], + false, + true // wildcards + ); for (int i=1; idict().found(entryNames[i])) { - entPtr = &entPtr->dict().lookupEntry(entryNames[i]); + entPtr = &entPtr->dict().lookupEntry + ( + entryNames[i], + false, + true // wildcards + ); } else { diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H index 201ce66811..71ac334133 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H @@ -151,18 +151,20 @@ public: const IOobject& fieldIoObject ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); }; diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C index 06f4f6a634..2ad0735ccf 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C @@ -131,7 +131,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -151,7 +151,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. label curF = cp[faceI] - 1; // Is the face on the boundary? @@ -282,7 +282,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField // It is necessary to create a copy of the addressing array to // take care of the face direction offset trick. - // + // { labelList curAddr(faceProcAddressing_[procI]); @@ -342,7 +342,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -452,11 +452,12 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField } -// Reconstruct and write all volume fields +// Reconstruct and write all/selected volume fields template void Foam::fvFieldReconstructor::reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -468,27 +469,29 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvVolumeField(*fieldIter())().write(); + reconstructFvVolumeField(*fieldIter())().write(); + } } - Info<< endl; } } -// Reconstruct and write all surface fields +// Reconstruct and write all/selected surface fields template void Foam::fvFieldReconstructor::reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -500,18 +503,19 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvSurfaceField(*fieldIter())().write(); + reconstructFvSurfaceField(*fieldIter())().write(); + } } - Info<< endl; } } diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C index db14efe3e5..50818f4db9 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C +++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C @@ -48,9 +48,17 @@ int main(int argc, char *argv[]) argList::noParallel(); timeSelector::addOptions(); # include "addRegionOption.H" + argList::validOptions.insert("fields", "\"(list of fields)\""); + # include "setRootCase.H" # include "createTime.H" + HashSet selectedFields; + if (args.options().found("fields")) + { + IStringStream(args.options()["fields"])() >> selectedFields; + } + // determine the processor count directly label nProcs = 0; while (dir(args.path()/(word("processor") + name(nProcs)))) @@ -184,13 +192,37 @@ int main(int argc, char *argv[]) procMeshes.boundaryProcAddressing() ); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); - fvReconstructor.reconstructFvSurfaceFields(objects); + fvReconstructor.reconstructFvSurfaceFields + ( + objects, + selectedFields + ); } else { diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamI.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamI.H index 8f9f682cdc..427c0a6fde 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamI.H +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamI.H @@ -37,7 +37,7 @@ inline Foam::word Foam::vtkPV3Foam::getFirstWord(const char* str) { ++n; } - return word(str, n); + return word(str, n, true); } else { diff --git a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C index b848e330ea..3f7a99040b 100644 --- a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C +++ b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C @@ -139,7 +139,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) IOobject LESPropertiesHeader ( - "RASProperties", + "LESProperties", runTime.constant(), mesh, IOobject::MUST_READ, diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict index 4c528ea243..e69c0386ef 100644 --- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict +++ b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict @@ -22,6 +22,7 @@ Ny 25 ); Nz 30; - + +symmetric true; // ************************************************************************* // diff --git a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C index f24729da97..ec5d019310 100644 --- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C +++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C @@ -164,7 +164,16 @@ int main(int argc, char *argv[]) forAll(dictList, i) { doneKeys[i] = dictList[i].keyword(); - dictList.set(i, fieldDict.lookupEntry(doneKeys[i]).clone()); + dictList.set + ( + i, + fieldDict.lookupEntry + ( + doneKeys[i], + false, + true + ).clone() + ); fieldDict.remove(doneKeys[i]); } // Add remaining entries diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C index aa30790b5c..f6312b32a4 100644 --- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C +++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C @@ -100,6 +100,11 @@ int main(int argc, char *argv[]) meshSubsetDict.lookup("addFaceNeighbours") ); + Switch invertSelection + ( + meshSubsetDict.lookup("invertSelection") + ); + // Mark the cells for the subset // Faces to subset @@ -337,6 +342,27 @@ int main(int argc, char *argv[]) << " faces because of addFaceNeighbours" << endl; } + + if (invertSelection) + { + Info<< "Inverting selection." << endl; + boolList newFacesToSubset(facesToSubset.size()); + + forAll(facesToSubset, i) + { + if (facesToSubset[i]) + { + newFacesToSubset[i] = false; + } + else + { + newFacesToSubset[i] = true; + } + } + facesToSubset.transfer(newFacesToSubset); + } + + // Create subsetted surface labelList pointMap; labelList faceMap; diff --git a/bin/mpirunDebug b/bin/mpirunDebug index 013a63997f..8d0be2375c 100755 --- a/bin/mpirunDebug +++ b/bin/mpirunDebug @@ -70,15 +70,17 @@ if [ ! -x "$exec" ]; then exit 1 fi +if [ ! "$PWD" ]; then + PWD=`pwd` +fi +echo "run $args" > $PWD/gdbCommands +echo "where" >> $PWD/gdbCommands +echo "Constructed gdb initialization file $PWD/gdbCommands" -echo "run $args" > $HOME/gdbCommands -echo "where" >> $HOME/gdbCommands -echo "Constructed gdb initialization file $HOME/gdbCommands" - -$ECHO "Choose running method: 1)gdb+xterm 2)gdb 3)log 4)log+xterm 5)xterm+valgrind: \c" +$ECHO "Choose running method: 0)normal 1)gdb+xterm 2)gdb 3)log 4)log+xterm 5)xterm+valgrind: \c" read method -if [ "$method" -ne 1 -a "$method" -ne 2 -a "$method" -ne 3 -a "$method" -ne 4 -a "$method" -ne 5 ]; then +if [ "$method" -ne 0 -a "$method" -ne 1 -a "$method" -ne 2 -a "$method" -ne 3 -a "$method" -ne 4 -a "$method" -ne 5 ]; then printUsage exit 1 fi @@ -119,15 +121,15 @@ fi echo "**sourceFoam:$sourceFoam" -rm -f $HOME/mpirun.schema -touch $HOME/mpirun.schema +rm -f $PWD/mpirun.schema +touch $PWD/mpirun.schema proc=0 xpos=0 ypos=0 for ((proc=0; proc<$nProcs; proc++)) do - procCmdFile="$HOME/processor${proc}.sh" + procCmdFile="$PWD/processor${proc}.sh" procLog="processor${proc}.log" geom="-geometry 120x20+$xpos+$ypos" node="" @@ -141,22 +143,25 @@ do fi echo "#!/bin/sh" > $procCmdFile - if [ "$method" -eq 1 ]; then - echo "$sourceFoam; cd $PWD; gdb -command $HOME/gdbCommands $exec 2>&1 | tee $procLog; read dummy" >> $procCmdFile + if [ "$method" -eq 0 ]; then + echo "$sourceFoam; cd $PWD; $exec $args | tee $procLog" >> $procCmdFile + echo "${node}$procCmdFile" >> $PWD/mpirun.schema + elif [ "$method" -eq 1 ]; then + echo "$sourceFoam; cd $PWD; gdb -command $PWD/gdbCommands $exec 2>&1 | tee $procLog; read dummy" >> $procCmdFile #echo "$sourceFoam; cd $PWD; $exec $args; read dummy" >> $procCmdFile - echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema + echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema elif [ "$method" -eq 2 ]; then - echo "$sourceFoam; cd $PWD; gdb -command $HOME/gdbCommands >& $procLog" >> $procCmdFile - echo "${node}$procCmdFile" >> $HOME/mpirun.schema + echo "$sourceFoam; cd $PWD; gdb -command $PWD/gdbCommands >& $procLog" >> $procCmdFile + echo "${node}$procCmdFile" >> $PWD/mpirun.schema elif [ "$method" -eq 3 ]; then echo "$sourceFoam; cd $PWD; $exec $args >& $procLog" >> $procCmdFile - echo "${node}$procCmdFile" >> $HOME/mpirun.schema + echo "${node}$procCmdFile" >> $PWD/mpirun.schema elif [ "$method" -eq 4 ]; then echo "$sourceFoam; cd $PWD; $exec $args 2>&1 | tee $procLog; read dummy" >> $procCmdFile - echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema + echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema elif [ "$method" -eq 5 ]; then echo "$sourceFoam; cd $PWD; valgrind $exec $args; read dummy" >> $procCmdFile - echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema + echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema fi chmod +x $procCmdFile @@ -176,10 +181,17 @@ do echo " tail -f $procLog" done -$ECHO "Constructed $HOME/mpirun.schema file. Press return to execute.\c" -read dummy +cmd="" if [ .$WM_MPLIB = .OPENMPI ]; then - mpirun -app $HOME/mpirun.schema & dfdy + scalarSquareMatrix& dfdy ) const = 0; }; diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C index fa00919a42..ec3db04574 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/KRR4/KRR4.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "KRR4.H" -#include "simpleMatrix.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -36,11 +35,11 @@ namespace Foam { addToRunTimeSelectionTable(ODESolver, KRR4, ODE); -const scalar +const scalar KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25, KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296; -const scalar +const scalar KRR4::gamma = 1.0/2.0, KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0, KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0, @@ -81,8 +80,8 @@ void Foam::KRR4::solve const ODE& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, + scalarField& dydx, + const scalar eps, const scalarField& yScale, const scalar hTry, scalar& hDid, @@ -109,14 +108,14 @@ void Foam::KRR4::solve a_[i][i] += 1.0/(gamma*h); } - simpleMatrix::LUDecompose(a_, pivotIndices_); + LUDecompose(a_, pivotIndices_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g1_); + LUBacksubstitute(a_, pivotIndices_, g1_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g2_); + LUBacksubstitute(a_, pivotIndices_, g2_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g3_); + LUBacksubstitute(a_, pivotIndices_, g3_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g4_); + LUBacksubstitute(a_, pivotIndices_, g4_); for (register label i=0; i dfdy_; - mutable Matrix a_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; mutable labelList pivotIndices_; static const int maxtry = 40; static const scalar safety, grow, pgrow, shrink, pshrink, errcon; - static const scalar + static const scalar gamma, a21, a31, a32, c21, c31, c32, c41, c42, c43, diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 75383bc38a..164b92d299 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.H @@ -60,8 +60,8 @@ class SIBS static const scalar safe1, safe2, redMax, redMin, scaleMX; mutable scalarField a_; - mutable Matrix alpha_; - mutable Matrix d_p_; + mutable scalarSquareMatrix alpha_; + mutable scalarSquareMatrix d_p_; mutable scalarField x_p_; mutable scalarField err_; @@ -69,7 +69,7 @@ class SIBS mutable scalarField ySeq_; mutable scalarField yErr_; mutable scalarField dfdx_; - mutable Matrix dfdy_; + mutable scalarSquareMatrix dfdy_; mutable label first_, kMax_, kOpt_; mutable scalar epsOld_, xNew_; @@ -82,7 +82,7 @@ void SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -97,7 +97,7 @@ void polyExtrapolate scalarField& yz, scalarField& dy, scalarField& x_p, - Matrix& d_p + scalarSquareMatrix& d_p ) const; diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C index f540f26d11..671f1acec9 100644 --- a/src/ODE/ODESolvers/SIBS/SIMPR.C +++ b/src/ODE/ODESolvers/SIBS/SIMPR.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "SIBS.H" -#include "simpleMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -36,7 +35,7 @@ void Foam::SIBS::SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -44,7 +43,7 @@ void Foam::SIBS::SIMPR { scalar h = deltaX/nSteps; - Matrix a(n_, n_); + scalarSquareMatrix a(n_); for (register label i=0; i::LUDecompose(a, pivotIndices); + LUDecompose(a, pivotIndices); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); scalarField del(yEnd); scalarField ytemp(n_); @@ -83,7 +82,7 @@ void Foam::SIBS::SIMPR yEnd[i] = h*yEnd[i] - del[i]; } - simpleMatrix::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i& d + scalarSquareMatrix& d ) const { label n = yz.size(); diff --git a/src/OSspecific/Unix/regularExpression.H b/src/OSspecific/Unix/regularExpression.H index dc574a4c9f..9924caef28 100644 --- a/src/OSspecific/Unix/regularExpression.H +++ b/src/OSspecific/Unix/regularExpression.H @@ -100,7 +100,7 @@ public: // Member functions //- Matches? - inline bool matches(const string& s) + inline bool matches(const string& s) const { size_t nmatch = 0; regmatch_t *pmatch = NULL; diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 7b137dedc0..e24fd6b5ec 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -39,6 +39,7 @@ $(strings)/word/word.C $(strings)/word/wordIO.C $(strings)/fileName/fileName.C $(strings)/fileName/fileNameIO.C +$(strings)/keyType/keyTypeIO.C primitives/random/Random.C @@ -177,8 +178,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C matrices/solution/solution.C -scalarMatrix = matrices/scalarMatrix -$(scalarMatrix)/scalarMatrix.C +scalarMatrices = matrices/scalarMatrices +$(scalarMatrices)/scalarMatrices.C +$(scalarMatrices)/SVD/SVD.C LUscalarMatrix = matrices/LUscalarMatrix $(LUscalarMatrix)/LUscalarMatrix.C diff --git a/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.C b/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.C index f0d3f8fd1a..577494aa1c 100644 --- a/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.C +++ b/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.C @@ -22,15 +22,31 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ +#include "word.H" #include "Ostream.H" #include "token.H" +#include "keyType.H" +#include "IOstreams.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +//- Write keyType +Foam::Ostream& Foam::Ostream::write(const keyType& s) +{ + // Write as word? + if (s.isWildCard()) + { + return write(static_cast(s)); + } + else + { + return write(static_cast(s)); + } +} + + //- Decrememt the indent level void Foam::Ostream::decrIndent() { @@ -47,7 +63,7 @@ void Foam::Ostream::decrIndent() // Write the keyword to the Ostream followed by appropriate indentation -Foam::Ostream& Foam::Ostream::writeKeyword(const Foam::word& keyword) +Foam::Ostream& Foam::Ostream::writeKeyword(const Foam::keyType& keyword) { indent(); write(keyword); diff --git a/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.H b/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.H index bcd081df56..6dc7df259a 100644 --- a/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.H +++ b/src/OpenFOAM/db/IOstreams/IOstreams/Ostream.H @@ -35,6 +35,7 @@ Description #define Ostream_H #include "IOstream.H" +#include "keyType.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -105,6 +106,9 @@ public: //- Write word virtual Ostream& write(const word&) = 0; + //- Write keyType + virtual Ostream& write(const keyType&); + //- Write string virtual Ostream& write(const string&) = 0; @@ -146,7 +150,7 @@ public: //- Write the keyword to the Ostream followed by // appropriate indentation - Ostream& writeKeyword(const word& keyword); + Ostream& writeKeyword(const keyType& keyword); // Stream state functions diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C index 6b078aa25d..7a1a680643 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C +++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "Pstream.H" diff --git a/src/OpenFOAM/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C index 8a38ff9ba2..ad960e7405 100644 --- a/src/OpenFOAM/db/dictionary/dictionary.C +++ b/src/OpenFOAM/db/dictionary/dictionary.C @@ -34,6 +34,72 @@ defineTypeNameAndDebug(Foam::dictionary, 0); const Foam::dictionary Foam::dictionary::null; + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::dictionary::findInWildcards +( + const bool wildCardMatch, + const word& Keyword, + DLList::const_iterator& wcLink, + DLList >::const_iterator& reLink +) const +{ + if (wildCardEntries_.size() > 0) + { + //wcLink = wildCardEntries_.begin(); + //reLink = wildCardRegexps_.end(); + + while (wcLink != wildCardEntries_.end()) + { + if (!wildCardMatch && wcLink()->keyword() == Keyword) + { + return true; + } + else if (wildCardMatch && reLink()->matches(Keyword)) + { + return true; + } + + ++reLink; + ++wcLink; + } + } + + return false; +} + + +bool Foam::dictionary::findInWildcards +( + const bool wildCardMatch, + const word& Keyword, + DLList::iterator& wcLink, + DLList >::iterator& reLink +) +{ + if (wildCardEntries_.size() > 0) + { + while (wcLink != wildCardEntries_.end()) + { + if (!wildCardMatch && wcLink()->keyword() == Keyword) + { + return true; + } + else if (wildCardMatch && reLink()->matches(Keyword)) + { + return true; + } + + ++reLink; + ++wcLink; + } + } + + return false; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::dictionary::dictionary() @@ -60,6 +126,18 @@ Foam::dictionary::dictionary ) { hashedEntries_.insert(iter().keyword(), &iter()); + + if (iter().keyword().isWildCard()) + { + wildCardEntries_.insert(&iter()); + wildCardRegexps_.insert + ( + autoPtr + ( + new regularExpression(iter().keyword()) + ) + ); + } } } @@ -81,6 +159,18 @@ Foam::dictionary::dictionary ) { hashedEntries_.insert(iter().keyword(), &iter()); + + if (iter().keyword().isWildCard()) + { + wildCardEntries_.insert(&iter()); + wildCardRegexps_.insert + ( + autoPtr + ( + new regularExpression(iter().keyword()) + ) + ); + } } } @@ -133,13 +223,29 @@ bool Foam::dictionary::found(const word& keyword, bool recursive) const { return true; } - else if (recursive && &parent_ != &dictionary::null) - { - return parent_.found(keyword, recursive); - } else { - return false; + if (wildCardEntries_.size() > 0) + { + DLList::const_iterator wcLink = wildCardEntries_.begin(); + DLList >::const_iterator reLink = + wildCardRegexps_.begin(); + + // Find in wildcards using regular expressions only + if (findInWildcards(true, keyword, wcLink, reLink)) + { + return true; + } + } + + if (recursive && &parent_ != &dictionary::null) + { + return parent_.found(keyword, recursive); + } + else + { + return false; + } } } @@ -147,16 +253,31 @@ bool Foam::dictionary::found(const word& keyword, bool recursive) const const Foam::entry* Foam::dictionary::lookupEntryPtr ( const word& keyword, - bool recursive + bool recursive, + bool wildCardMatch ) const { HashTable::const_iterator iter = hashedEntries_.find(keyword); if (iter == hashedEntries_.end()) { + if (wildCardMatch && wildCardEntries_.size() > 0) + { + DLList::const_iterator wcLink = + wildCardEntries_.begin(); + DLList >::const_iterator reLink = + wildCardRegexps_.begin(); + + // Find in wildcards using regular expressions only + if (findInWildcards(wildCardMatch, keyword, wcLink, reLink)) + { + return wcLink(); + } + } + if (recursive && &parent_ != &dictionary::null) { - return parent_.lookupEntryPtr(keyword, recursive); + return parent_.lookupEntryPtr(keyword, recursive, wildCardMatch); } else { @@ -171,19 +292,34 @@ const Foam::entry* Foam::dictionary::lookupEntryPtr Foam::entry* Foam::dictionary::lookupEntryPtr ( const word& keyword, - bool recursive + bool recursive, + bool wildCardMatch ) { HashTable::iterator iter = hashedEntries_.find(keyword); if (iter == hashedEntries_.end()) { + if (wildCardMatch && wildCardEntries_.size() > 0) + { + DLList::iterator wcLink = + wildCardEntries_.begin(); + DLList >::iterator reLink = + wildCardRegexps_.begin(); + // Find in wildcards using regular expressions only + if (findInWildcards(wildCardMatch, keyword, wcLink, reLink)) + { + return wcLink(); + } + } + if (recursive && &parent_ != &dictionary::null) { return const_cast(parent_).lookupEntryPtr ( keyword, - recursive + recursive, + wildCardMatch ); } else @@ -199,16 +335,17 @@ Foam::entry* Foam::dictionary::lookupEntryPtr const Foam::entry& Foam::dictionary::lookupEntry ( const word& keyword, - bool recursive + bool recursive, + bool wildCardMatch ) const { - const entry* entryPtr = lookupEntryPtr(keyword, recursive); + const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch); if (entryPtr == NULL) { FatalIOErrorIn ( - "dictionary::lookupEntry(const word& keyword) const", + "dictionary::lookupEntry(const word&, bool, bool) const", *this ) << "keyword " << keyword << " is undefined in dictionary " << name() @@ -222,16 +359,18 @@ const Foam::entry& Foam::dictionary::lookupEntry Foam::ITstream& Foam::dictionary::lookup ( const word& keyword, - bool recursive + bool recursive, + bool wildCardMatch ) const { - return lookupEntry(keyword, recursive).stream(); + return lookupEntry(keyword, recursive, wildCardMatch).stream(); } bool Foam::dictionary::isDict(const word& keyword) const { - const entry* entryPtr = lookupEntryPtr(keyword); + // Find non-recursive with wildcards + const entry* entryPtr = lookupEntryPtr(keyword, false, true); if (entryPtr) { @@ -246,7 +385,7 @@ bool Foam::dictionary::isDict(const word& keyword) const const Foam::dictionary* Foam::dictionary::subDictPtr(const word& keyword) const { - const entry* entryPtr = lookupEntryPtr(keyword); + const entry* entryPtr = lookupEntryPtr(keyword, false, true); if (entryPtr) { @@ -261,7 +400,8 @@ const Foam::dictionary* Foam::dictionary::subDictPtr(const word& keyword) const const Foam::dictionary& Foam::dictionary::subDict(const word& keyword) const { - const entry* entryPtr = lookupEntryPtr(keyword); + const entry* entryPtr = lookupEntryPtr(keyword, false, true); + if (entryPtr == NULL) { FatalIOErrorIn @@ -278,7 +418,8 @@ const Foam::dictionary& Foam::dictionary::subDict(const word& keyword) const Foam::dictionary& Foam::dictionary::subDict(const word& keyword) { - entry* entryPtr = lookupEntryPtr(keyword); + entry* entryPtr = lookupEntryPtr(keyword, false, true); + if (entryPtr == NULL) { FatalIOErrorIn @@ -314,7 +455,10 @@ Foam::wordList Foam::dictionary::toc() const bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry) { - HashTable::iterator iter = hashedEntries_.find(entryPtr->keyword()); + HashTable::iterator iter = hashedEntries_.find + ( + entryPtr->keyword() + ); if (mergeEntry && iter != hashedEntries_.end()) { @@ -336,6 +480,19 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry) if (hashedEntries_.insert(entryPtr->keyword(), entryPtr)) { entryPtr->name() = name_ + "::" + entryPtr->keyword(); + + if (entryPtr->keyword().isWildCard()) + { + wildCardEntries_.insert(entryPtr); + wildCardRegexps_.insert + ( + autoPtr + ( + new regularExpression(entryPtr->keyword()) + ) + ); + } + return true; } else @@ -356,6 +513,18 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry) entryPtr->name() = name_ + "::" + entryPtr->keyword(); IDLList::append(entryPtr); + if (entryPtr->keyword().isWildCard()) + { + wildCardEntries_.insert(entryPtr); + wildCardRegexps_.insert + ( + autoPtr + ( + new regularExpression(entryPtr->keyword()) + ) + ); + } + return true; } else @@ -376,27 +545,37 @@ void Foam::dictionary::add(const entry& e, bool mergeEntry) add(e.clone(*this).ptr(), mergeEntry); } -void Foam::dictionary::add(const word& k, const word& w, bool overwrite) +void Foam::dictionary::add(const keyType& k, const word& w, bool overwrite) { add(new primitiveEntry(k, token(w)), overwrite); } -void Foam::dictionary::add(const word& k, const Foam::string& s, bool overwrite) +void Foam::dictionary::add +( + const keyType& k, + const Foam::string& s, + bool overwrite +) { add(new primitiveEntry(k, token(s)), overwrite); } -void Foam::dictionary::add(const word& k, const label l, bool overwrite) +void Foam::dictionary::add(const keyType& k, const label l, bool overwrite) { add(new primitiveEntry(k, token(l)), overwrite); } -void Foam::dictionary::add(const word& k, const scalar s, bool overwrite) +void Foam::dictionary::add(const keyType& k, const scalar s, bool overwrite) { add(new primitiveEntry(k, token(s)), overwrite); } -void Foam::dictionary::add(const word& k, const dictionary& d, bool mergeEntry) +void Foam::dictionary::add +( + const keyType& k, + const dictionary& d, + bool mergeEntry +) { add(new dictionaryEntry(k, *this, d), mergeEntry); } @@ -404,7 +583,7 @@ void Foam::dictionary::add(const word& k, const dictionary& d, bool mergeEntry) void Foam::dictionary::set(entry* entryPtr) { - entry* existingPtr = lookupEntryPtr(entryPtr->keyword()); + entry* existingPtr = lookupEntryPtr(entryPtr->keyword(), false, true); // clear dictionary so merge acts like overwrite if (existingPtr && existingPtr->isDict()) @@ -420,7 +599,7 @@ void Foam::dictionary::set(const entry& e) set(e.clone(*this).ptr()); } -void Foam::dictionary::set(const word& k, const dictionary& d) +void Foam::dictionary::set(const keyType& k, const dictionary& d) { set(new dictionaryEntry(k, *this, d)); } @@ -432,6 +611,19 @@ bool Foam::dictionary::remove(const word& Keyword) if (iter != hashedEntries_.end()) { + // Delete from wildcards first + DLList::iterator wcLink = + wildCardEntries_.begin(); + DLList >::iterator reLink = + wildCardRegexps_.begin(); + + // Find in wildcards using exact match only + if (findInWildcards(false, Keyword, wcLink, reLink)) + { + wildCardEntries_.remove(wcLink); + wildCardRegexps_.remove(reLink); + } + IDLList::remove(iter()); delete iter(); hashedEntries_.erase(iter); @@ -447,8 +639,8 @@ bool Foam::dictionary::remove(const word& Keyword) bool Foam::dictionary::changeKeyword ( - const word& oldKeyword, - const word& newKeyword, + const keyType& oldKeyword, + const keyType& newKeyword, bool forceOverwrite ) { @@ -466,6 +658,18 @@ bool Foam::dictionary::changeKeyword return false; } + if (iter()->keyword().isWildCard()) + { + FatalErrorIn + ( + "dictionary::changeKeyword(const word&, const word&, bool)" + ) << "Old keyword "<< oldKeyword + << " is a wildcard." + << "Wildcard replacement not yet implemented." + << exit(FatalError); + } + + HashTable::iterator iter2 = hashedEntries_.find(newKeyword); // newKeyword already exists @@ -473,14 +677,33 @@ bool Foam::dictionary::changeKeyword { if (forceOverwrite) { + if (iter2()->keyword().isWildCard()) + { + // Delete from wildcards first + DLList::iterator wcLink = + wildCardEntries_.begin(); + DLList >::iterator reLink = + wildCardRegexps_.begin(); + + // Find in wildcards using exact match only + if (findInWildcards(false, iter2()->keyword(), wcLink, reLink)) + { + wildCardEntries_.remove(wcLink); + wildCardRegexps_.remove(reLink); + } + } + IDLList::replace(iter2(), iter()); delete iter2(); hashedEntries_.erase(iter2); + } else { - WarningIn("dictionary::changeKeyword(const word&, const word&)") - << "cannot rename keyword "<< oldKeyword + WarningIn + ( + "dictionary::changeKeyword(const word&, const word&, bool)" + ) << "cannot rename keyword "<< oldKeyword << " to existing keyword " << newKeyword << " in dictionary " << name() << endl; return false; @@ -493,6 +716,18 @@ bool Foam::dictionary::changeKeyword hashedEntries_.erase(oldKeyword); hashedEntries_.insert(newKeyword, iter()); + if (newKeyword.isWildCard()) + { + wildCardEntries_.insert(iter()); + wildCardRegexps_.insert + ( + autoPtr + ( + new regularExpression(newKeyword) + ) + ); + } + return true; } @@ -579,6 +814,7 @@ void Foam::dictionary::operator=(const dictionary& dict) // Create clones of the entries in the given dictionary // resetting the parentDict to this dictionary + for ( IDLList::const_iterator iter = dict.begin(); @@ -586,17 +822,7 @@ void Foam::dictionary::operator=(const dictionary& dict) ++iter ) { - IDLList::append(iter().clone(*this).ptr()); - } - - for - ( - IDLList::iterator iter = begin(); - iter != end(); - ++iter - ) - { - hashedEntries_.insert(iter().keyword(), &iter()); + add(iter().clone(*this).ptr()); } } diff --git a/src/OpenFOAM/db/dictionary/dictionary.H b/src/OpenFOAM/db/dictionary/dictionary.H index 05872cce61..ac6ca65389 100644 --- a/src/OpenFOAM/db/dictionary/dictionary.H +++ b/src/OpenFOAM/db/dictionary/dictionary.H @@ -27,7 +27,12 @@ Class Description A list of keyword definitions, which are a keyword followed by any number - of values (e.g. words and numbers). + of values (e.g. words and numbers). The keywords can represent wildcards + which are matched using Posix regular expressions. The general order for + searching is + - exact match + - wildcard match (in reverse order) + - optional recursion into subdictionaries The dictionary class is the base class for IOdictionary. It also serves as a bootstrap dictionary for the objectRegistry data @@ -49,11 +54,13 @@ SourceFiles #include "entry.H" #include "IDLList.H" +#include "DLList.H" #include "fileName.H" #include "ITstream.H" #include "HashTable.H" #include "wordList.H" #include "className.H" +#include "regularExpression.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -80,12 +87,40 @@ class dictionary //- Dictionary name fileName name_; - //- HashTable of the enries held on the DL-list for quick lookup + //- HashTable of the entries held on the DL-list for quick lookup HashTable hashedEntries_; //- Parent dictionary const dictionary& parent_; + //- Wildcard entries + DLList wildCardEntries_; + + //- Wildcard precompiled regex + DLList > wildCardRegexps_; + + // Private Member Functions + + //- Search wildcard table either for exact match or for regular + // expression match. + bool findInWildcards + ( + const bool wildCardMatch, + const word& Keyword, + DLList::const_iterator& wcLink, + DLList >::const_iterator& reLink + ) const; + + //- Search wildcard table either for exact match or for regular + // expression match. + bool findInWildcards + ( + const bool wildCardMatch, + const word& Keyword, + DLList::iterator& wcLink, + DLList >::iterator& reLink + ); + public: @@ -181,24 +216,44 @@ public: //- Find and return an entry data stream pointer if present // otherwise return NULL. - // If recursive search parent dictionaries + // If recursive search parent dictionaries. If wildCardMatch + // use wildcards. const entry* lookupEntryPtr ( - const word&, bool recursive=false + const word&, + bool recursive, + bool wildCardMatch ) const; //- Find and return an entry data stream pointer for manipulation // if present otherwise return NULL. - // If recursive search parent dictionaries - entry* lookupEntryPtr(const word&, bool recursive=false); + // If recursive search parent dictionaries. If wildCardMatch + // use wildcards. + entry* lookupEntryPtr + ( + const word&, + bool recursive, + bool wildCardMatch + ); //- Find and return an entry data stream if present otherwise error. - // If recursive search parent dictionaries - const entry& lookupEntry(const word&, bool recursive=false) const; + // If recursive search parent dictionaries. If wildCardMatch + // use wildcards. + const entry& lookupEntry + ( + const word&, + bool recursive, + bool wildCardMatch + ) const; //- Find and return an entry data stream // If recursive search parent dictionaries - ITstream& lookup(const word&, bool recursive=false) const; + ITstream& lookup + ( + const word&, + bool recursive=false, + bool wildCardMatch=true + ) const; //- Find and return a T, // if not found return the given default value @@ -208,7 +263,8 @@ public: ( const word&, const T&, - bool recursive=false + bool recursive=false, + bool wildCardMatch=true ) const; //- Find and return a T, if not found return the given @@ -219,7 +275,8 @@ public: ( const word&, const T&, - bool recursive=false + bool recursive=false, + bool wildCardMatch=true ); //- Find an entry if present, and assign to T @@ -229,7 +286,8 @@ public: ( const word&, T&, - bool recursive=false + bool recursive=false, + bool wildCardMatch=true ) const; //- Check if entry is a sub-dictionary @@ -248,7 +306,6 @@ public: //- Return the table of contents wordList toc() const; - // Editing //- Add a new entry @@ -263,25 +320,25 @@ public: //- Add a word entry // optionally overwrite an existing entry - void add(const word& keyword, const word&, bool overwrite=false); + void add(const keyType&, const word&, bool overwrite=false); //- Add a string entry // optionally overwrite an existing entry - void add(const word& keyword, const string&, bool overwrite=false); + void add(const keyType&, const string&, bool overwrite=false); //- Add a label entry // optionally overwrite an existing entry - void add(const word& keyword, const label, bool overwrite=false); + void add(const keyType&, const label, bool overwrite=false); //- Add a scalar entry // optionally overwrite an existing entry - void add (const word& keyword, const scalar, bool overwrite=false); + void add (const keyType&, const scalar, bool overwrite=false); //- Add a dictionary entry // optionally merge with an existing sub-dictionary void add ( - const word& keyword, + const keyType& keyword, const dictionary&, bool mergeEntry=false ); @@ -289,7 +346,7 @@ public: //- Add a T entry // optionally overwrite an existing entry template - void add(const word& keyword, const T&, bool overwrite=false); + void add(const keyType& keyword, const T&, bool overwrite=false); //- Assign a new entry, overwrite any existing entry void set(entry*); @@ -298,11 +355,11 @@ public: void set(const entry&); //- Assign a dictionary entry, overwrite any existing entry - void set(const word& keyword, const dictionary&); + void set(const keyType& keyword, const dictionary&); //- Assign a T entry, overwrite any existing entry template - void set(const word& keyword, const T&); + void set(const keyType& keyword, const T&); //- Remove an entry specified by keyword bool remove(const word& keyword); @@ -311,8 +368,8 @@ public: // optionally forcing overwrite of an existing entry bool changeKeyword ( - const word& oldKeyword, - const word& newKeyword, + const keyType& oldKeyword, + const keyType& newKeyword, bool forceOverwrite = false ); @@ -361,11 +418,13 @@ public: // Global Operators -//- Combine dictionaries starting from the entries in dict1 and then including those from dict2. +//- Combine dictionaries starting from the entries in dict1 and then including +// those from dict2. // Warn, but do not overwrite the entries from dict1. dictionary operator+(const dictionary& dict1, const dictionary& dict2); -//- Combine dictionaries starting from the entries in dict1 and then including those from dict2. +//- Combine dictionaries starting from the entries in dict1 and then including +// those from dict2. // Do not overwrite the entries from dict1. dictionary operator|(const dictionary& dict1, const dictionary& dict2); diff --git a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C index 9beb835bcb..6dea212507 100644 --- a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C +++ b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C @@ -30,7 +30,7 @@ License Foam::dictionaryEntry::dictionaryEntry ( - const word& key, + const keyType& key, const dictionary& parentDict, const dictionary& dict ) diff --git a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.H b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.H index 2812c1b0bb..1909f4851d 100644 --- a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.H +++ b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.H @@ -79,7 +79,7 @@ public: //- Construct from the keyword, parent dictionary and a Istream dictionaryEntry ( - const word& keyword, + const keyType& keyword, const dictionary& parentDict, Istream& is ); @@ -87,7 +87,7 @@ public: //- Construct from the keyword, parent dictionary and a dictionary dictionaryEntry ( - const word& keyword, + const keyType& keyword, const dictionary& parentDict, const dictionary& dict ); diff --git a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntryIO.C b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntryIO.C index b6c2c2ebce..9cce3eb7de 100644 --- a/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntryIO.C +++ b/src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntryIO.C @@ -27,7 +27,9 @@ Description \*---------------------------------------------------------------------------*/ +#include "keyType.H" #include "dictionaryEntry.H" +#include "IOstreams.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -43,14 +45,14 @@ Foam::dictionaryEntry::dictionaryEntry is.fatalCheck ( "dictionaryEntry::dictionaryEntry" - "(Istream& is, const dictionary& parentDict)" + "(const dictionary& parentDict, Istream& is)" ); } Foam::dictionaryEntry::dictionaryEntry ( - const word& key, + const keyType& key, const dictionary& parentDict, Istream& is ) @@ -63,7 +65,7 @@ Foam::dictionaryEntry::dictionaryEntry is.fatalCheck ( "dictionaryEntry::dictionaryEntry" - "(const word& keyword, const dictionary& parentDict, Istream& is)" + "(const keyType& keyword, const dictionary& parentDict, Istream& is)" ); } diff --git a/src/OpenFOAM/db/dictionary/dictionaryIO.C b/src/OpenFOAM/db/dictionary/dictionaryIO.C index 53069af0b0..b1d5fbcbc4 100644 --- a/src/OpenFOAM/db/dictionary/dictionaryIO.C +++ b/src/OpenFOAM/db/dictionary/dictionaryIO.C @@ -71,7 +71,7 @@ bool Foam::dictionary::substituteKeyword(const word& keyword) word varName = keyword(1, keyword.size()-1); // lookup the variable name in the given dictionary.... - const entry* ePtr = lookupEntryPtr(varName, true); + const entry* ePtr = lookupEntryPtr(varName, true, true); // ...if defined insert its entries into this dictionary... if (ePtr != NULL) @@ -137,6 +137,8 @@ Foam::Istream& Foam::operator>>(Istream& is, dictionary& dict) dict.clear(); dict.hashedEntries_.clear(); + dict.wildCardEntries_.clear(); + dict.wildCardRegexps_.clear(); dict.read(is); return is; diff --git a/src/OpenFOAM/db/dictionary/dictionaryTemplates.C b/src/OpenFOAM/db/dictionary/dictionaryTemplates.C index 3dc32a7b51..a47a976814 100644 --- a/src/OpenFOAM/db/dictionary/dictionaryTemplates.C +++ b/src/OpenFOAM/db/dictionary/dictionaryTemplates.C @@ -34,10 +34,11 @@ T Foam::dictionary::lookupOrDefault ( const word& keyword, const T& deflt, - bool recursive + bool recursive, + bool wildCardMatch ) const { - const entry* entryPtr = lookupEntryPtr(keyword, recursive); + const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch); if (entryPtr == NULL) { @@ -55,10 +56,11 @@ T Foam::dictionary::lookupOrAddDefault ( const word& keyword, const T& deflt, - bool recursive + bool recursive, + bool wildCardMatch ) { - const entry* entryPtr = lookupEntryPtr(keyword, recursive); + const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch); if (entryPtr == NULL) { @@ -77,10 +79,11 @@ bool Foam::dictionary::readIfPresent ( const word& k, T& val, - bool recursive + bool recursive, + bool wildCardMatch ) const { - const entry* entryPtr = lookupEntryPtr(k, recursive); + const entry* entryPtr = lookupEntryPtr(k, recursive, wildCardMatch); if (entryPtr == NULL) { @@ -95,16 +98,17 @@ bool Foam::dictionary::readIfPresent template -void Foam::dictionary::add(const word& k, const T& t, bool overwrite) +void Foam::dictionary::add(const keyType& k, const T& t, bool overwrite) { add(new primitiveEntry(k, t), overwrite); } template -void Foam::dictionary::set(const word& k, const T& t) +void Foam::dictionary::set(const keyType& k, const T& t) { set(new primitiveEntry(k, t)); } + // ************************************************************************* // diff --git a/src/OpenFOAM/db/dictionary/entry/entry.C b/src/OpenFOAM/db/dictionary/entry/entry.C index 93cf677b0b..54581d869b 100644 --- a/src/OpenFOAM/db/dictionary/entry/entry.C +++ b/src/OpenFOAM/db/dictionary/entry/entry.C @@ -30,7 +30,7 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::entry::entry(const word& keyword) +Foam::entry::entry(const keyType& keyword) : keyword_(keyword) {} diff --git a/src/OpenFOAM/db/dictionary/entry/entry.H b/src/OpenFOAM/db/dictionary/entry/entry.H index 5afab95b5a..5ed8b929b6 100644 --- a/src/OpenFOAM/db/dictionary/entry/entry.H +++ b/src/OpenFOAM/db/dictionary/entry/entry.H @@ -42,6 +42,7 @@ SourceFiles #ifndef entry_H #define entry_H +#include "keyType.H" #include "IDLList.H" #include "fileName.H" #include "autoPtr.H" @@ -70,13 +71,13 @@ class entry // Private data //- Keyword of entry - word keyword_; + keyType keyword_; // Private Member Functions //- Get the next valid keyword otherwise return false - static bool getKeyword(word& keyword, Istream& is); + static bool getKeyword(keyType& keyword, Istream& is); public: @@ -84,7 +85,7 @@ public: // Constructors //- Construct from keyword - entry(const word& keyword); + entry(const keyType& keyword); //- Construct as copy entry(const entry&); @@ -116,13 +117,13 @@ public: // Member functions //- Return keyword - const word& keyword() const + const keyType& keyword() const { return keyword_; } //- Return non-const access to keyword - word& keyword() + keyType& keyword() { return keyword_; } diff --git a/src/OpenFOAM/db/dictionary/entry/entryIO.C b/src/OpenFOAM/db/dictionary/entry/entryIO.C index 830c21cf8a..0142b302ad 100644 --- a/src/OpenFOAM/db/dictionary/entry/entryIO.C +++ b/src/OpenFOAM/db/dictionary/entry/entryIO.C @@ -32,7 +32,7 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -bool Foam::entry::getKeyword(word& keyword, Istream& is) +bool Foam::entry::getKeyword(keyType& keyword, Istream& is) { token keywordToken; @@ -57,6 +57,12 @@ bool Foam::entry::getKeyword(word& keyword, Istream& is) keyword = keywordToken.wordToken(); return true; } + else if (keywordToken.isString()) + { + // Enable wildcards + keyword = keywordToken.stringToken(); + return true; + } // If it is the end of the dictionary or file return false... else if (keywordToken == token::END_BLOCK || is.eof()) { @@ -67,7 +73,7 @@ bool Foam::entry::getKeyword(word& keyword, Istream& is) { cerr<< "--> FOAM Warning : " << std::endl << " From function " - << "entry::getKeyword(word& keyword, Istream& is)" << std::endl + << "entry::getKeyword(keyType& keyword, Istream& is)" << std::endl << " in file " << __FILE__ << " at line " << __LINE__ << std::endl << " Reading " << is.name().c_str() << std::endl @@ -84,7 +90,7 @@ bool Foam::entry::New(dictionary& parentDict, Istream& is) { is.fatalCheck("entry::New(const dictionary& parentDict, Istream& is)"); - word keyword; + keyType keyword; // Get the next keyword and if invalid return false if (!getKeyword(keyword, is)) @@ -115,7 +121,13 @@ bool Foam::entry::New(dictionary& parentDict, Istream& is) // Deal with duplicate entries bool mergeEntry = false; - entry* existingPtr = parentDict.lookupEntryPtr(keyword); + // See (using exact match) if entry already present + entry* existingPtr = parentDict.lookupEntryPtr + ( + keyword, + false, + false + ); if (existingPtr) { if (functionEntries::inputModeEntry::overwrite()) @@ -158,7 +170,7 @@ Foam::autoPtr Foam::entry::New(Istream& is) { is.fatalCheck("entry::New(Istream& is)"); - word keyword; + keyType keyword; // Get the next keyword and if invalid return false if (!getKeyword(keyword, is)) diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C index df888cd064..bc77cefef3 100644 --- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C +++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C @@ -29,7 +29,7 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::primitiveEntry::primitiveEntry(const word& key, const ITstream& tokens) +Foam::primitiveEntry::primitiveEntry(const keyType& key, const ITstream& tokens) : entry(key), ITstream(tokens) @@ -38,7 +38,7 @@ Foam::primitiveEntry::primitiveEntry(const word& key, const ITstream& tokens) } -Foam::primitiveEntry::primitiveEntry(const word& keyword, const token& t) +Foam::primitiveEntry::primitiveEntry(const keyType& keyword, const token& t) : entry(keyword), ITstream(keyword, tokenList(1, t)) @@ -47,7 +47,7 @@ Foam::primitiveEntry::primitiveEntry(const word& keyword, const token& t) Foam::primitiveEntry::primitiveEntry ( - const word& keyword, + const keyType& keyword, const tokenList& tokens ) : diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.H b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.H index b97452818a..86d8afd61d 100644 --- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.H +++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.H @@ -108,23 +108,23 @@ public: // Constructors //- Construct from keyword and a Istream - primitiveEntry(const word& keyword, Istream&); + primitiveEntry(const keyType& keyword, Istream&); //- Construct from keyword, parent dictionary and a Istream - primitiveEntry(const word& keyword, const dictionary&, Istream&); + primitiveEntry(const keyType& keyword, const dictionary&, Istream&); //- Construct from keyword and a ITstream - primitiveEntry(const word& keyword, const ITstream&); + primitiveEntry(const keyType& keyword, const ITstream&); //- Construct from keyword and a token - primitiveEntry(const word&, const token&); + primitiveEntry(const keyType&, const token&); //- Construct from keyword and a tokenList - primitiveEntry(const word&, const tokenList&); + primitiveEntry(const keyType&, const tokenList&); //- Construct from keyword and a T template - primitiveEntry(const word&, const T&); + primitiveEntry(const keyType&, const T&); autoPtr clone(const dictionary&) const { diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C index 10bc3c78c8..58cf475303 100644 --- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C +++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C @@ -81,7 +81,7 @@ bool Foam::primitiveEntry::expandVariable word varName = w(1, w.size()-1); // lookup the variable name in the given dictionary.... - const entry* ePtr = dict.lookupEntryPtr(varName, true); + const entry* ePtr = dict.lookupEntryPtr(varName, true, true); // ...if defined insert its tokens into this if (ePtr != NULL) @@ -218,7 +218,7 @@ void Foam::primitiveEntry::readEntry(const dictionary& dict, Istream& is) Foam::primitiveEntry::primitiveEntry ( - const word& key, + const keyType& key, const dictionary& dict, Istream& is ) @@ -236,7 +236,7 @@ Foam::primitiveEntry::primitiveEntry } -Foam::primitiveEntry::primitiveEntry(const word& key, Istream& is) +Foam::primitiveEntry::primitiveEntry(const keyType& key, Istream& is) : entry(key), ITstream diff --git a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryTemplates.C b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryTemplates.C index fce303173b..38e35b054a 100644 --- a/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryTemplates.C +++ b/src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryTemplates.C @@ -30,7 +30,7 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::primitiveEntry::primitiveEntry(const word& keyword, const T& t) +Foam::primitiveEntry::primitiveEntry(const keyType& keyword, const T& t) : entry(keyword), ITstream(keyword, tokenList(10)) diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C new file mode 100644 index 0000000000..a31058eace --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -0,0 +1,100 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 "DiagonalMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +: + List(min(a.n(), a.m())) +{ + forAll(*this, i) + { + this->operator[](i) = a[i][i]; + } +} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size) +: + List(size) +{} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size, const Type& val) +: + List(size, val) +{} + + +template +Foam::DiagonalMatrix& Foam::DiagonalMatrix::invert() +{ + forAll(*this, i) + { + Type x = this->operator[](i); + if (mag(x) < VSMALL) + { + this->operator[](i) = Type(0); + } + else + { + this->operator[](i) = Type(1)/x; + } + } + + return this; +} + + +template +Foam::DiagonalMatrix Foam::inv(const DiagonalMatrix& A) +{ + DiagonalMatrix Ainv = A; + + forAll(A, i) + { + Type x = A[i]; + if (mag(x) < VSMALL) + { + Ainv[i] = Type(0); + } + else + { + Ainv[i] = Type(1)/x; + } + } + + return Ainv; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H new file mode 100644 index 0000000000..ae8ceffc1c --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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::DiagonalMatrix + +Description + DiagonalMatrix is a 2D diagonal matrix of objects + of type Type, size nxn + +SourceFiles + DiagonalMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef DiagonalMatrix_H +#define DiagonalMatrix_H + +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // + +template class Matrix; + +/*---------------------------------------------------------------------------*\ + Class DiagonalMatrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class DiagonalMatrix +: + public List +{ +public: + + // Constructors + + //- Construct from diagonal component of a Matrix + template + DiagonalMatrix(const Matrix&); + + //- Construct empty from size + DiagonalMatrix(const label size); + + //- Construct from size and a value + DiagonalMatrix(const label, const Type&); + + + // Member functions + + //- Invert the diaganol matrix and return itself + DiagonalMatrix& invert(); +}; + + +// Global functions + +//- Return the diagonal Matrix inverse +template +DiagonalMatrix inv(const DiagonalMatrix&); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "DiagonalMatrix.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 5de8693f2f..3e8fc3fe0f 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -31,9 +31,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::LUscalarMatrix::LUscalarMatrix(const Matrix& matrix) +Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), pivotIndices_(n()) { LUDecompose(*this, pivotIndices_); @@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix nCells += lduMatrices[i].size(); } - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(lduMatrices); } @@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix else { label nCells = ldum.lduAddr().size(); - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(ldum, interfaceCoeffs, interfaces); } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index d5e79faf19..2255cee236 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef LUscalarMatrix_H #define LUscalarMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "labelList.H" #include "FieldField.H" #include "lduInterfaceFieldPtrsList.H" @@ -55,7 +55,7 @@ class procLduMatrix; class LUscalarMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data @@ -78,7 +78,7 @@ class LUscalarMatrix void convert(const PtrList& lduMatrices); - //- Print the ratio of the mag-sum of the off-diagonal coefficients + //- Print the ratio of the mag-sum of the off-diagonal coefficients // to the mag-diagonal void printDiagonalDominance() const; @@ -87,8 +87,8 @@ public: // Constructors - //- Construct from Matrix and perform LU decomposition - LUscalarMatrix(const Matrix&); + //- Construct from scalarSquareMatrix and perform LU decomposition + LUscalarMatrix(const scalarSquareMatrix&); //- Construct from lduMatrix and perform LU decomposition LUscalarMatrix diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C index cbb59e2c86..ac0bc17207 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C @@ -28,16 +28,16 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::LUscalarMatrix::solve(Field& sourceSol) const +template +void Foam::LUscalarMatrix::solve(Field& sourceSol) const { if (Pstream::parRun()) { - Field completeSourceSol(n()); + Field completeSourceSol(n()); if (Pstream::master()) { - typename Field::subField + typename Field::subField ( completeSourceSol, sourceSol.size() @@ -58,7 +58,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } @@ -77,7 +77,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const { LUBacksubstitute(*this, pivotIndices_, completeSourceSol); - sourceSol = typename Field::subField + sourceSol = typename Field::subField ( completeSourceSol, sourceSol.size() @@ -98,7 +98,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } diff --git a/src/OpenFOAM/containers/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C similarity index 59% rename from src/OpenFOAM/containers/Matrix/Matrix.C rename to src/OpenFOAM/matrices/Matrix/Matrix.C index 733e04034e..577cbcabcb 100644 --- a/src/OpenFOAM/containers/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -28,13 +28,13 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -template -void Foam::Matrix::allocate() +template +void Foam::Matrix::allocate() { if (n_ && m_) { - v_ = new T*[n_]; - v_[0] = new T[n_*m_]; + v_ = new Type*[n_]; + v_[0] = new Type[n_*m_]; for (register label i=1; i::allocate() // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // -template -Foam::Matrix::~Matrix() +template +Foam::Matrix::~Matrix() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } } @@ -59,16 +59,16 @@ Foam::Matrix::~Matrix() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -const Foam::Matrix& Foam::Matrix::null() +template +const Foam::Matrix& Foam::Matrix::null() { - Matrix* nullPtr = reinterpret_cast*>(NULL); + Matrix* nullPtr = reinterpret_cast*>(NULL); return *nullPtr; } -template -Foam::Matrix::Matrix(const label n, const label m) +template +Foam::Matrix::Matrix(const label n, const label m) : n_(n), m_(m), @@ -76,7 +76,7 @@ Foam::Matrix::Matrix(const label n, const label m) { if (n_ < 0 || m_ < 0) { - FatalErrorIn("Matrix::Matrix(const label n, const label m)") + FatalErrorIn("Matrix::Matrix(const label n, const label m)") << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -85,8 +85,8 @@ Foam::Matrix::Matrix(const label n, const label m) } -template -Foam::Matrix::Matrix(const label n, const label m, const T& a) +template +Foam::Matrix::Matrix(const label n, const label m, const Type& a) : n_(n), m_(m), @@ -96,7 +96,7 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) { FatalErrorIn ( - "Matrix::Matrix(const label n, const label m, const T&)" + "Matrix::Matrix(const label n, const label m, const T&)" ) << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -105,7 +105,7 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) if (v_) { - T* v = v_[0]; + Type* v = v_[0]; label nm = n_*m_; @@ -117,8 +117,8 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) } -template -Foam::Matrix::Matrix(const Matrix& a) +template +Foam::Matrix::Matrix(const Matrix& a) : n_(a.n_), m_(a.m_), @@ -127,8 +127,8 @@ Foam::Matrix::Matrix(const Matrix& a) if (a.v_) { allocate(); - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::Matrix(const Matrix& a) } } - -template -void Foam::Matrix::clear() + +template +void Foam::Matrix::clear() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } n_ = 0; @@ -153,8 +153,8 @@ void Foam::Matrix::clear() } -template -void Foam::Matrix::transfer(Matrix& a) +template +void Foam::Matrix::transfer(Matrix& a) { clear(); @@ -169,14 +169,32 @@ void Foam::Matrix::transfer(Matrix& a) } +template +Form Foam::Matrix::T() const +{ + const Matrix& A = *this; + Form At(m(), n()); + + for (register label i=0; i -void Foam::Matrix::operator=(const T& t) +template +void Foam::Matrix::operator=(const Type& t) { if (v_) { - T* v = v_[0]; + Type* v = v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const T& t) // Assignment operator. Takes linear time. -template -void Foam::Matrix::operator=(const Matrix& a) +template +void Foam::Matrix::operator=(const Matrix& a) { if (this == &a) { - FatalErrorIn("Matrix::operator=(const Matrix&)") + FatalErrorIn("Matrix::operator=(const Matrix&)") << "attempted assignment to self" << abort(FatalError); } @@ -204,12 +222,12 @@ void Foam::Matrix::operator=(const Matrix& a) n_ = a.n_; m_ = a.m_; allocate(); - } + } if (v_) { - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const Matrix& a) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // -template -const T& Foam::max(const Matrix& a) +template +const Type& Foam::max(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMaxI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("max(const Matrix&)") + FatalErrorIn("max(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -254,15 +272,15 @@ const T& Foam::max(const Matrix& a) } -template -const T& Foam::min(const Matrix& a) +template +const Type& Foam::min(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMinI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("min(const Matrix&)") + FatalErrorIn("min(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -288,15 +306,15 @@ const T& Foam::min(const Matrix& a) // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // -template -Foam::Matrix Foam::operator-(const Matrix& a) +template +Form Foam::operator-(const Matrix& a) { - Matrix na(a.n_, a.m_); + Form na(a.n_, a.m_); if (a.n_ && a.m_) { - T* nav = na.v_[0]; - const T* av = a.v_[0]; + Type* nav = na.v_[0]; + const Type* av = a.v_[0]; label nm = a.n_*a.m_; for (register label i=0; i Foam::operator-(const Matrix& a) } -template -Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) +template +Form Foam::operator+(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator+(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) +template +Form Foam::operator-(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator-(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator*(const scalar s, const Matrix& a) +template +Form Foam::operator*(const scalar s, const Matrix& a) { - Matrix sa(a.n_, a.m_); + Form sa(a.n_, a.m_); if (a.n_ && a.m_) { - T* sav = sa.v_[0]; - const T* av = a.v_[0]; + Type* sav = sa.v_[0]; + const Type* av = a.v_[0]; - label nm = a.n_*a.m_;; + label nm = a.n_*a.m_; for (register label i=0; i class Matrix; +template class Matrix; -template const T& max(const Matrix&); -template const T& min(const Matrix&); +template Istream& operator>> +( + Istream&, + Matrix& +); -template Matrix operator-(const Matrix&); -template Matrix operator+(const Matrix&, const Matrix&); -template Matrix operator-(const Matrix&, const Matrix&); -template Matrix operator*(const scalar, const Matrix&); - -template Istream& operator>>(Istream&, Matrix&); -template Ostream& operator<<(Ostream&, const Matrix&); +template Ostream& operator<< +( + Ostream&, + const Matrix& +); /*---------------------------------------------------------------------------*\ Class Matrix Declaration \*---------------------------------------------------------------------------*/ -template +template class Matrix { // Private data @@ -78,7 +79,7 @@ class Matrix label n_, m_; //- Row pointers - T** __restrict__ v_; + Type** __restrict__ v_; //- Allocate the storage for the row-pointers and the data // and set the row pointers @@ -97,16 +98,16 @@ public: //- Construct with given number of rows and columns // and value for all elements. - Matrix(const label n, const label m, const T&); + Matrix(const label n, const label m, const Type&); //- Copy constructor. - Matrix(const Matrix&); + Matrix(const Matrix&); //- Construct from Istream. Matrix(Istream&); //- Clone - inline autoPtr > clone() const; + inline autoPtr > clone() const; // Destructor @@ -117,7 +118,7 @@ public: // Member functions //- Return a null Matrix - static const Matrix& null(); + static const Matrix& null(); // Access @@ -148,48 +149,64 @@ public: //- Transfer the contents of the argument Matrix into this Matrix // and annull the argument Matrix. - void transfer(Matrix&); + void transfer(Matrix&); + + + //- Return the transpose of the matrix + Form T() const; // Member operators //- Return subscript-checked element of Matrix. - inline T* operator[](const label); + inline Type* operator[](const label); //- Return subscript-checked element of constant Matrix. - inline const T* operator[](const label) const; + inline const Type* operator[](const label) const; //- Assignment operator. Takes linear time. - void operator=(const Matrix&); + void operator=(const Matrix&); //- Assignment of all entries to the given value - void operator=(const T&); - - - // Friend functions - - friend const T& max(const Matrix&); - friend const T& min(const Matrix&); - - - // Friend operators - - friend Matrix operator-(const Matrix&); - friend Matrix operator+(const Matrix&, const Matrix&); - friend Matrix operator-(const Matrix&, const Matrix&); - friend Matrix operator*(const scalar, const Matrix&); + void operator=(const Type&); // IOstream operators //- Read Matrix from Istream, discarding contents of existing Matrix. - friend Istream& operator>> (Istream&, Matrix&); + friend Istream& operator>> (Istream&, Matrix&); // Write Matrix to Ostream. - friend Ostream& operator<< (Ostream&, const Matrix&); + friend Ostream& operator<< (Ostream&, const Matrix&); }; +// Global functions and operators + +template const Type& max(const Matrix&); +template const Type& min(const Matrix&); + +template Form operator-(const Matrix&); + +template Form operator+ +( + const Matrix&, + const Matrix& +); + +template Form operator- +( + const Matrix&, + const Matrix& +); + +template Form operator* +( + const scalar, + const Matrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/containers/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H similarity index 66% rename from src/OpenFOAM/containers/Matrix/MatrixI.H rename to src/OpenFOAM/matrices/Matrix/MatrixI.H index d0e3fe8fdf..8d758cf5b1 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -26,8 +26,8 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -inline Foam::Matrix::Matrix() +template +inline Foam::Matrix::Matrix() : n_(0), m_(0), @@ -35,71 +35,67 @@ inline Foam::Matrix::Matrix() {} -template -inline Foam::autoPtr > Foam::Matrix::clone() const +template +inline Foam::autoPtr > Foam::Matrix::clone() const { - return autoPtr >(new Matrix(*this)); + return autoPtr >(new Matrix(*this)); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // //- Return the number of rows -template -inline Foam::label Foam::Matrix::n() const +template +inline Foam::label Foam::Matrix::n() const { return n_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::m() const +template +inline Foam::label Foam::Matrix::m() const { return m_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::size() const +template +inline Foam::label Foam::Matrix::size() const { return n_*m_; } -// Check index i is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checki(const label i) const +template +inline void Foam::Matrix::checki(const label i) const { if (!n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "attempt to access element from zero sized row" << abort(FatalError); } else if (i<0 || i>=n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "index " << i << " out of range 0 ... " << n_-1 << abort(FatalError); } } -// Check index j is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checkj(const label j) const +template +inline void Foam::Matrix::checkj(const label j) const { if (!m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "attempt to access element from zero sized column" << abort(FatalError); } else if (j<0 || j>=m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "index " << j << " out of range 0 ... " << m_-1 << abort(FatalError); } @@ -108,9 +104,8 @@ inline void Foam::Matrix::checkj(const label j) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -// Return subscript-checked element access -template -inline T* Foam::Matrix::operator[](const label i) +template +inline Type* Foam::Matrix::operator[](const label i) { # ifdef FULLDEBUG checki(i); @@ -119,9 +114,8 @@ inline T* Foam::Matrix::operator[](const label i) } -// Return subscript-checked const element access -template -inline const T* Foam::Matrix::operator[](const label i) const +template +inline const Type* Foam::Matrix::operator[](const label i) const { # ifdef FULLDEBUG checki(i); diff --git a/src/OpenFOAM/containers/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C similarity index 83% rename from src/OpenFOAM/containers/Matrix/MatrixIO.C rename to src/OpenFOAM/matrices/Matrix/MatrixIO.C index 710e42876e..15f0d5345e 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixIO.C +++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C @@ -32,8 +32,8 @@ License // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // -template -Foam::Matrix::Matrix(Istream& is) +template +Foam::Matrix::Matrix(Istream& is) : n_(0), m_(0), @@ -43,17 +43,17 @@ Foam::Matrix::Matrix(Istream& is) } -template -Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) +template +Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) { // Anull matrix M.clear(); - is.fatalCheck("operator>>(Istream&, Matrix&)"); + is.fatalCheck("operator>>(Istream&, Matrix&)"); token firstToken(is); - is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); + is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); if (firstToken.isLabel()) { @@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) label nm = M.n_*M.m_; // Read list contents depending on data format - if (is.format() == IOstream::ASCII || !contiguous()) + if (is.format() == IOstream::ASCII || !contiguous()) { // Read beginning of contents char listDelimiter = is.readBeginList("Matrix"); @@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) if (nm) { M.allocate(); - T* v = M.v_[0]; + Type* v = M.v_[0]; if (listDelimiter == token::BEGIN_LIST) { @@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading entry" ); } @@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - T element; + Type element; is >> element; is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the single entry" ); @@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) if (nm) { M.allocate(); - T* v = M.v_[0]; + Type* v = M.v_[0]; - is.read(reinterpret_cast(v), nm*sizeof(T)); + is.read(reinterpret_cast(v), nm*sizeof(Type)); is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the binary block" ); } @@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) + FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) << "incorrect first token, expected , found " << firstToken.info() << exit(FatalIOError); @@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } -template -Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { label nm = M.n_*M.m_; os << M.n() << token::SPACE << M.m(); // Write list contents depending on data format - if (os.format() == IOstream::ASCII || !contiguous()) + if (os.format() == IOstream::ASCII || !contiguous()) { if (nm) { bool uniform = false; - const T* v = M.v_[0]; + const Type* v = M.v_[0]; - if (nm > 1 && contiguous()) + if (nm > 1 && contiguous()) { uniform = true; @@ -187,7 +187,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) // Write end of contents delimiter os << token::END_BLOCK; } - else if (nm < 10 && contiguous()) + else if (nm < 10 && contiguous()) { // Write size of list and start contents delimiter os << token::BEGIN_LIST; @@ -246,7 +246,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { if (nm) { - os.write(reinterpret_cast(M.v_[0]), nm*sizeof(T)); + os.write(reinterpret_cast(M.v_[0]), nm*sizeof(Type)); } } diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H new file mode 100644 index 0000000000..e5e6b16777 --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H @@ -0,0 +1,92 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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::RectangularMatrix + +Description + A templated 2D rectangular matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. + +SourceFiles + RectangularMatrixI.H + RectangularMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RectangularMatrix_H +#define RectangularMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class RectangularMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline RectangularMatrix(); + + //- Construct given number of rows and columns, + inline RectangularMatrix(const label m, const label n); + + //- Construct with given number of rows and columns + // and value for all elements. + inline RectangularMatrix(const label m, const label n, const Type&); + + //- Construct from Istream. + inline RectangularMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "RectangularMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H new file mode 100644 index 0000000000..1cacd2f11a --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::RectangularMatrix::RectangularMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n +) +: + Matrix, Type>(m, n) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n, + const Type& t +) +: + Matrix, Type>(m, n, t) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::RectangularMatrix::clone() const +{ + return autoPtr >(new RectangularMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H new file mode 100644 index 0000000000..27160dc239 --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -0,0 +1,97 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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::SquareMatrix + +Description + A templated 2D square matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. + +SourceFiles + SquareMatrixI.H + SquareMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SquareMatrix_H +#define SquareMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class SquareMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline SquareMatrix(); + + //- Construct given number of rows/columns. + inline SquareMatrix(const label n); + + //- Construct given number of rows and columns, + // It checks that m == n. + inline SquareMatrix(const label m, const label n); + + //- Construct with given number of rows and rows + // and value for all elements. + // It checks that m == n. + inline SquareMatrix(const label m, const label n, const Type&); + + //- Construct from Istream. + inline SquareMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "SquareMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H new file mode 100644 index 0000000000..d3ee1cd6bf --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::SquareMatrix::SquareMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label n) +: + Matrix, Type>(n, n) +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label m, const label n) +: + Matrix, Type>(m, n) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix(const label m, const label n)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} + +template +inline Foam::SquareMatrix::SquareMatrix +( + const label m, + const label n, + const Type& t +) +: + Matrix, Type>(m, n, t) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix" + "(const label m, const label n, const Type&)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} + +template +inline Foam::SquareMatrix::SquareMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::SquareMatrix::clone() const +{ + return autoPtr >(new SquareMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C new file mode 100644 index 0000000000..13ef29ed9b --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C @@ -0,0 +1,403 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2005 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 "SVD.H" +#include "scalarList.H" +#include "scalarMatrices.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * // + +Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) +: + U_(A), + V_(A.m(), A.m()), + S_(A.m()), + VSinvUt_(A.m(), A.n()), + nZeros_(0) +{ + // SVDcomp to find U_, V_ and S_ - the singular values + + const label Um = U_.m(); + const label Un = U_.n(); + + scalarList rv1(Um); + scalar g = 0; + scalar scale = 0; + scalar s = 0; + scalar anorm = 0; + label l = 0; + + for (label i = 0; i < Um; i++) + { + l = i+2; + rv1[i] = scale*g; + g = s = scale = 0; + + if (i < Un) + { + for (label k = i; k < Un; k++) + { + scale += mag(U_[k][i]); + } + + if (scale != 0) + { + for (label k = i; k < Un; k++) + { + U_[k][i] /= scale; + s += U_[k][i]*U_[k][i]; + } + + scalar f = U_[i][i]; + g = -sign(Foam::sqrt(s), f); + scalar h = f*g - s; + U_[i][i] = f - g; + + for (label j = l-1; j < Um; j++) + { + s = 0; + for (label k = i; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + f = s/h; + for (label k = i; k < A.n(); k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label k = i; k < Un; k++) + { + U_[k][i] *= scale; + } + } + } + + S_[i] = scale*g; + + g = s = scale = 0; + + if (i+1 <= Un && i != Um) + { + for (label k = l-1; k < Um; k++) + { + scale += mag(U_[i][k]); + } + + if (scale != 0) + { + for (label k=l-1; k < Um; k++) + { + U_[i][k] /= scale; + s += U_[i][k]*U_[i][k]; + } + + scalar f = U_[i][l-1]; + g = -sign(Foam::sqrt(s),f); + scalar h = f*g - s; + U_[i][l-1] = f - g; + + for (label k = l-1; k < Um; k++) + { + rv1[k] = U_[i][k]/h; + } + + for (label j = l-1; j < Un; j++) + { + s = 0; + for (label k = l-1; k < Um; k++) + { + s += U_[j][k]*U_[i][k]; + } + + for (label k = l-1; k < Um; k++) + { + U_[j][k] += s*rv1[k]; + } + } + for (label k = l-1; k < Um; k++) + { + U_[i][k] *= scale; + } + } + } + + anorm = max(anorm, mag(S_[i]) + mag(rv1[i])); + } + + for (label i = Um-1; i >= 0; i--) + { + if (i < Um-1) + { + if (g != 0) + { + for (label j = l; j < Um; j++) + { + V_[j][i] = (U_[i][j]/U_[i][l])/g; + } + + for (label j=l; j < Um; j++) + { + s = 0; + for (label k = l; k < Um; k++) + { + s += U_[i][k]*V_[k][j]; + } + + for (label k = l; k < Um; k++) + { + V_[k][j] += s*V_[k][i]; + } + } + } + + for (label j = l; j < Um;j++) + { + V_[i][j] = V_[j][i] = 0.0; + } + } + + V_[i][i] = 1; + g = rv1[i]; + l = i; + } + + for (label i = min(Um, Un) - 1; i >= 0; i--) + { + l = i+1; + g = S_[i]; + + for (label j = l; j < Um; j++) + { + U_[i][j] = 0.0; + } + + if (g != 0) + { + g = 1.0/g; + + for (label j = l; j < Um; j++) + { + s = 0; + for (label k = l; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + scalar f = (s/U_[i][i])*g; + + for (label k = i; k < Un; k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label j = i; j < Un; j++) + { + U_[j][i] *= g; + } + } + else + { + for (label j = i; j < Un; j++) + { + U_[j][i] = 0.0; + } + } + + ++U_[i][i]; + } + + for (label k = Um-1; k >= 0; k--) + { + for (label its = 0; its < 35; its++) + { + bool flag = true; + + label nm; + for (l = k; l >= 0; l--) + { + nm = l-1; + if (mag(rv1[l]) + anorm == anorm) + { + flag = false; + break; + } + if (mag(S_[nm]) + anorm == anorm) break; + } + + if (flag) + { + scalar c = 0.0; + s = 1.0; + for (label i = l-1; i < k+1; i++) + { + scalar f = s*rv1[i]; + rv1[i] = c*rv1[i]; + + if (mag(f) + anorm == anorm) break; + + g = S_[i]; + scalar h = sqrtSumSqr(f, g); + S_[i] = h; + h = 1.0/h; + c = g*h; + s = -f*h; + + for (label j = 0; j < Un; j++) + { + scalar y = U_[j][nm]; + scalar z = U_[j][i]; + U_[j][nm] = y*c + z*s; + U_[j][i] = z*c - y*s; + } + } + } + + scalar z = S_[k]; + + if (l == k) + { + if (z < 0.0) + { + S_[k] = -z; + + for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k]; + } + break; + } + if (its == 34) + { + WarningIn + ( + "SVD::SVD" + "(scalarRectangularMatrix& A, const scalar minCondition)" + ) << "no convergence in 35 SVD iterations" + << endl; + } + + scalar x = S_[l]; + nm = k-1; + scalar y = S_[nm]; + g = rv1[nm]; + scalar h = rv1[k]; + scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y); + g = sqrtSumSqr(f, 1.0); + f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x; + scalar c = 1.0; + s = 1.0; + + for (label j = l; j <= nm; j++) + { + label i = j + 1; + g = rv1[i]; + y = S_[i]; + h = s*g; + g = c*g; + scalar z = sqrtSumSqr(f, h); + rv1[j] = z; + c = f/z; + s = h/z; + f = x*c + g*s; + g = g*c - x*s; + h = y*s; + y *= c; + + for (label jj = 0; jj < Um; jj++) + { + x = V_[jj][j]; + z = V_[jj][i]; + V_[jj][j] = x*c + z*s; + V_[jj][i] = z*c - x*s; + } + + z = sqrtSumSqr(f, h); + S_[j] = z; + if (z) + { + z = 1.0/z; + c = f*z; + s = h*z; + } + f = c*g + s*y; + x = c*y - s*g; + + for (label jj=0; jj < Un; jj++) + { + y = U_[jj][j]; + z = U_[jj][i]; + U_[jj][j] = y*c + z*s; + U_[jj][i] = z*c - y*s; + } + } + rv1[l] = 0.0; + rv1[k] = f; + S_[k] = x; + } + } + + // zero singular values that are less than minCondition*maxS + const scalar minS = minCondition*S_[findMax(S_)]; + for (label i = 0; i < S_.size(); i++) + { + if (S_[i] <= minS) + { + //Info << "Removing " << S_[i] << " < " << minS << endl; + S_[i] = 0; + nZeros_++; + } + } + + // now multiply out to find the pseudo inverse of A, VSinvUt_ + multiply(VSinvUt_, V_, inv(S_), U_.T()); + + // test SVD + /*scalarRectangularMatrix SVDA(A.n(), A.m()); + multiply(SVDA, U_, S_, transpose(V_)); + scalar maxDiff = 0; + scalar diff = 0; + for(label i = 0; i < A.n(); i++) + { + for(label j = 0; j < A.m(); j++) + { + diff = mag(A[i][j] - SVDA[i][j]); + if (diff > maxDiff) maxDiff = diff; + } + } + Info << "Maximum discrepancy between A and svd(A) = " << maxDiff << endl; + + if (maxDiff > 4) + { + Info << "singular values " << S_ << endl; + } + */ +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H new file mode 100644 index 0000000000..48e581198a --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2005 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::SVD + +Description + Singular value decomposition of a rectangular matrix. + +SourceFiles + SVDI.H + SVD.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SVD_H +#define SVD_H + +#include "scalarMatrices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +/*---------------------------------------------------------------------------*\ + Class SVD Declaration +\*---------------------------------------------------------------------------*/ + +class SVD +{ + // Private data + + //- Rectangular matrix with the same dimensions as the input + scalarRectangularMatrix U_; + + //- square matrix V + scalarRectangularMatrix V_; + + //- The singular values + DiagonalMatrix S_; + + //- The matrix product V S^(-1) U^T + scalarRectangularMatrix VSinvUt_; + + //- The number of zero singular values + label nZeros_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + SVD(const SVD&); + + //- Disallow default bitwise assignment + void operator=(const SVD&); + + template + inline const T sign(const T& a, const T& b); + + +public: + + // Constructors + + //- Construct from a rectangular Matrix + SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0); + + + // Access functions + + //- Return U + inline const scalarRectangularMatrix& U() const; + + //- Return the square matrix V + inline const scalarRectangularMatrix& V() const; + + //- Return the singular values + inline const scalarDiagonalMatrix& S() const; + + //- Return VSinvUt (the pseudo inverse) + inline const scalarRectangularMatrix& VSinvUt() const; + + //- Return the number of zero singular values + inline label nZeros() const; + + //- Return the minimum non-zero singular value + inline scalar minNonZeroS() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "SVDI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H new file mode 100644 index 0000000000..298aac1aeb --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // + +template +inline const T Foam::SVD::sign(const T& a, const T& b) +{ + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const +{ + return U_; +} + +inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const +{ + return V_; +} + +inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const +{ + return S_; +} + +inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const +{ + return VSinvUt_; +} + +inline Foam::label Foam::SVD::nZeros() const +{ + return nZeros_; +} + +inline Foam::scalar Foam::SVD::minNonZeroS() const +{ + scalar minS = S_[0]; + for(label i = 1; i < S_.size(); i++) + { + scalar s = S_[i]; + if (s > VSMALL && s < minS) minS = s; + } + return minS; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C new file mode 100644 index 0000000000..8ddd01a82c --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -0,0 +1,293 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 "scalarMatrices.H" +#include "SVD.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices +) +{ + label n = matrix.n(); + scalar vv[n]; + + for (register label i=0; i largestCoeff) + { + largestCoeff = temp; + } + } + + if (largestCoeff == 0.0) + { + FatalErrorIn + ( + "LUdecompose" + "(scalarSquareMatrix& matrix, labelList& rowIndices)" + ) << "Singular matrix" << exit(FatalError); + } + + vv[i] = 1.0/largestCoeff; + } + + for (register label j=0; j= largestCoeff) + { + largestCoeff = temp; + iMax = i; + } + } + + pivotIndices[j] = iMax; + + if (j != iMax) + { + scalar* __restrict__ matrixiMax = matrix[iMax]; + + for (register label k=0; k& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" + ) << "A and B must have identical inner dimensions but A.m = " + << A.m() << " and B.n = " << B.size() + << abort(FatalError); + } + + if (B.size() != C.n()) + { + FatalErrorIn + ( + "multiply(" + "const scalarRectangularMatrix& A, " + "const DiagonalMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" + ) << "B and C must have identical inner dimensions but B.m = " + << B.size() << " and C.n = " << C.n() + << abort(FatalError); + } + + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); + + for(register label i = 0; i < A.n(); i++) + { + for(register label g = 0; g < C.m(); g++) + { + for(register label l = 0; l < C.n(); l++) + { + ans[i][g] += C[l][g] * A[i][l]*B[l]; + } + } + } +} + + +Foam::RectangularMatrix Foam::SVDinv +( + const scalarRectangularMatrix& A, + scalar minCondition +) +{ + SVD svd(A, minCondition); + return svd.VSinvUt(); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H new file mode 100644 index 0000000000..17d8727516 --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 + scalarMatrices + +Description + Scalar matrices + +SourceFiles + scalarMatrices.C + scalarMatricesTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarMatrices_H +#define scalarMatrices_H + +#include "RectangularMatrix.H" +#include "SquareMatrix.H" +#include "DiagonalMatrix.H" +#include "scalarField.H" +#include "labelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +typedef RectangularMatrix scalarRectangularMatrix; +typedef SquareMatrix scalarSquareMatrix; +typedef DiagonalMatrix scalarDiagonalMatrix; + +//- Solve the matrix using Gaussian elimination with pivoting, +// returning the solution in the source +template +void solve(scalarSquareMatrix& matrix, Field& source); + +//- Solve the matrix using Gaussian elimination with pivoting +// and return the solution +template +void solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +); + +//- LU decompose the matrix with pivoting +void LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices +); + +//- LU back-substitution with given source, returning the solution +// in the source +template +void LUBacksubstitute +( + const scalarSquareMatrix& luMmatrix, + const labelList& pivotIndices, + Field& source +); + +//- Solve the matrix using LU decomposition with pivoting +// returning the LU form of the matrix and the solution in the source +template +void LUsolve(scalarSquareMatrix& matrix, Field& source); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B, + const scalarRectangularMatrix& C +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const DiagonalMatrix& B, + const scalarRectangularMatrix& C +); + +//- Return the inverse of matrix A using SVD +scalarRectangularMatrix SVDinv +( + const scalarRectangularMatrix& A, + scalar minCondition = 0 +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "scalarMatricesTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C similarity index 83% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C index b729e32fcc..131856b3fc 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C @@ -24,16 +24,16 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "Swap.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::scalarMatrix::solve +template +void Foam::solve ( - Matrix& tmpMatrix, - Field& sourceSol + scalarSquareMatrix& tmpMatrix, + Field& sourceSol ) { label n = tmpMatrix.n(); @@ -68,7 +68,7 @@ void Foam::scalarMatrix::solve // Check that the system of equations isn't singular if (mag(tmpMatrix[i][i]) < 1e-20) { - FatalErrorIn("scalarMatrix::solve()") + FatalErrorIn("solve(scalarSquareMatrix&, Field& sourceSol)") << "Singular Matrix" << exit(FatalError); } @@ -89,7 +89,7 @@ void Foam::scalarMatrix::solve // Back-substitution for (register label j=n-1; j>=0; j--) { - T ntempvec = pTraits::zero; + Type ntempvec = pTraits::zero; for (register label k=j+1; k -void Foam::scalarMatrix::solve(Field& psi, const Field& source) const +template +void Foam::solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +) { - Matrix tmpMatrix = *this; + scalarSquareMatrix tmpMatrix = matrix; psi = source; solve(tmpMatrix, psi); } -template -void Foam::scalarMatrix::LUBacksubstitute +template +void Foam::LUBacksubstitute ( - const Matrix& luMatrix, + const scalarSquareMatrix& luMatrix, const labelList& pivotIndices, - Field& sourceSol + Field& sourceSol ) { label n = luMatrix.n(); @@ -125,7 +130,7 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=0; i::zero) + else if (sum != pTraits::zero) { ii = i+1; } @@ -146,11 +151,11 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=n-1; i>=0; i--) { - T sum = sourceSol[i]; + Type sum = sourceSol[i]; const scalar* __restrict__ luMatrixi = luMatrix[i]; for (register label j=i+1; j -void Foam::scalarMatrix::LUsolve +template +void Foam::LUsolve ( - Matrix& matrix, - Field& sourceSol + scalarSquareMatrix& matrix, + Field& sourceSol ) { labelList pivotIndices(matrix.n()); diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C deleted file mode 100644 index 2ff9e39f18..0000000000 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C +++ /dev/null @@ -1,161 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 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 "scalarMatrix.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::scalarMatrix::scalarMatrix() -{} - - -Foam::scalarMatrix::scalarMatrix(const label mSize) -: - Matrix(mSize, mSize, 0.0) -{} - - -Foam::scalarMatrix::scalarMatrix(const Matrix& matrix) -: - Matrix(matrix) -{} - - -Foam::scalarMatrix::scalarMatrix(Istream& is) -: - Matrix(is) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::scalarMatrix::LUDecompose -( - Matrix& matrix, - labelList& pivotIndices -) -{ - label n = matrix.n(); - scalar vv[n]; - - for (register label i=0; i largestCoeff) - { - largestCoeff = temp; - } - } - - if (largestCoeff == 0.0) - { - FatalErrorIn - ( - "scalarMatrix::LUdecompose" - "(Matrix& matrix, labelList& rowIndices)" - ) << "Singular matrix" << exit(FatalError); - } - - vv[i] = 1.0/largestCoeff; - } - - for (register label j=0; j= largestCoeff) - { - largestCoeff = temp; - iMax = i; - } - } - - pivotIndices[j] = iMax; - - if (j != iMax) - { - scalar* __restrict__ matrixiMax = matrix[iMax]; - - for (register label k=0; k -Foam::simpleMatrix::simpleMatrix(const label mSize) +template +Foam::simpleMatrix::simpleMatrix(const label mSize) : - scalarMatrix(mSize), - source_(mSize, pTraits::zero) + scalarSquareMatrix(mSize), + source_(mSize, pTraits::zero) {} -template -Foam::simpleMatrix::simpleMatrix +template +Foam::simpleMatrix::simpleMatrix ( - const scalarMatrix& matrix, - const Field& source + const scalarSquareMatrix& matrix, + const Field& source ) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), source_(source) {} -template -Foam::simpleMatrix::simpleMatrix(Istream& is) +template +Foam::simpleMatrix::simpleMatrix(Istream& is) : - scalarMatrix(is), + scalarSquareMatrix(is), source_(is) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -Foam::Field Foam::simpleMatrix::solve() const +template +Foam::Field Foam::simpleMatrix::solve() const { - scalarMatrix tmpMatrix = *this; - Field sourceSol = source_; + scalarSquareMatrix tmpMatrix = *this; + Field sourceSol = source_; - scalarMatrix::solve(tmpMatrix, sourceSol); + Foam::solve(tmpMatrix, sourceSol); return sourceSol; } -template -Foam::Field Foam::simpleMatrix::LUsolve() const +template +Foam::Field Foam::simpleMatrix::LUsolve() const { - scalarMatrix luMatrix = *this; - Field sourceSol = source_; + scalarSquareMatrix luMatrix = *this; + Field sourceSol = source_; - scalarMatrix::LUsolve(luMatrix, sourceSol); + Foam::LUsolve(luMatrix, sourceSol); return sourceSol; } @@ -84,82 +84,82 @@ Foam::Field Foam::simpleMatrix::LUsolve() const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -void Foam::simpleMatrix::operator=(const simpleMatrix& m) +template +void Foam::simpleMatrix::operator=(const simpleMatrix& m) { if (this == &m) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Attempted assignment to self" << abort(FatalError); } if (n() != m.n()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size matrices" << abort(FatalError); } if (source_.size() != m.source_.size()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size source vectors" << abort(FatalError); } - scalarMatrix::operator=(m); + scalarSquareMatrix::operator=(m); source_ = m.source_; } // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // -template -Foam::simpleMatrix Foam::operator+ +template +Foam::simpleMatrix Foam::operator+ ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( - static_cast(m1) - + static_cast(m2), + static_cast(m1) + + static_cast(m2), m1.source_ + m2.source_ ); } -template -Foam::simpleMatrix Foam::operator- +template +Foam::simpleMatrix Foam::operator- ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( - static_cast(m1) - - static_cast(m2), + static_cast(m1) + - static_cast(m2), m1.source_ - m2.source_ ); } -template -Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) +template +Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) { - return simpleMatrix(s*m.matrix_, s*m.source_); + return simpleMatrix(s*m.matrix_, s*m.source_); } // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -template -Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) { - os << static_cast(m) << nl << m.source_; + os << static_cast(m) << nl << m.source_; return os; } diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H index d1c411452b..a972c787c8 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef simpleMatrix_H #define simpleMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -45,35 +45,14 @@ namespace Foam // Forward declaration of friend functions and operators -template +template class simpleMatrix; -template -simpleMatrix operator+ -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator- -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator* -( - const scalar, - const simpleMatrix& -); - -template +template Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); @@ -81,14 +60,14 @@ Ostream& operator<< Class simpleMatrix Declaration \*---------------------------------------------------------------------------*/ -template +template class simpleMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data - Field source_; + Field source_; public: @@ -99,25 +78,25 @@ public: simpleMatrix(const label); //- Construct from components - simpleMatrix(const scalarMatrix&, const Field&); + simpleMatrix(const scalarSquareMatrix&, const Field&); //- Construct from Istream simpleMatrix(Istream&); //- Construct as copy - simpleMatrix(const simpleMatrix&); + simpleMatrix(const simpleMatrix&); // Member Functions // Access - Field& source() + Field& source() { return source_; } - const Field& source() const + const Field& source() const { return source_; } @@ -125,49 +104,52 @@ public: //- Solve the matrix using Gaussian elimination with pivoting // and return the solution - Field solve() const; + Field solve() const; //- Solve the matrix using LU decomposition with pivoting // and return the solution - Field LUsolve() const; + Field LUsolve() const; // Member Operators - void operator=(const simpleMatrix&); - - - // Friend Operators - - friend simpleMatrix operator+ - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator- - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator* - ( - const scalar, - const simpleMatrix& - ); + void operator=(const simpleMatrix&); // Ostream Operator - friend Ostream& operator<< + friend Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); }; +// Global operators + +template +simpleMatrix operator+ +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator- +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator* +( + const scalar, + const simpleMatrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/primitives/strings/keyType/keyType.H similarity index 52% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H rename to src/OpenFOAM/primitives/strings/keyType/keyType.H index f2b0764f1a..4d4c358d5d 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H +++ b/src/OpenFOAM/primitives/strings/keyType/keyType.H @@ -23,88 +23,107 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::scalarMatrix + Foam::keyType Description - Foam::scalarMatrix + A class for handling keywords in dictionaries. + + A keyType is the keyword of a dictionary. It differs from word in that + it accepts wildcards. SourceFiles - scalarMatrix.C + keyType.C + keyTypeIO.C \*---------------------------------------------------------------------------*/ -#ifndef scalarMatrix_H -#define scalarMatrix_H +#ifndef keyType_H +#define keyType_H -#include "Matrix.H" -#include "scalarField.H" -#include "labelList.H" +#include "word.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { +// Forward declaration of classes +class Istream; +class Ostream; + + /*---------------------------------------------------------------------------*\ - Class scalarMatrix Declaration + Class keyType Declaration \*---------------------------------------------------------------------------*/ -class scalarMatrix +class keyType : - public Matrix + public word { + // Private member data + + bool isWildCard_; + + // Private Member Functions + + //- Disallow assignments where we cannot determine string/word type + void operator=(const std::string&); public: + // Constructors //- Construct null - scalarMatrix(); + inline keyType(); - //- Construct given size - scalarMatrix(const label); + //- Construct as copy + inline keyType(const keyType& s); - //- Construct from Matrix - scalarMatrix(const Matrix&); + //- Construct as copy of word + inline keyType(const word& s); + + //- Construct as copy of string. Expect it to be regular expression. + inline keyType(const string& s); + + //- Construct as copy of character array + inline keyType(const char* s); + + //- Construct as copy of std::string + inline keyType(const std::string& s, const bool isWildCard); //- Construct from Istream - scalarMatrix(Istream&); + keyType(Istream& is); - // Member Functions + // Member functions - //- Solve the matrix using Gaussian elimination with pivoting, - // returning the solution in the source - template - static void solve(Matrix& matrix, Field& source); + //- Is this character valid for a keyType + inline static bool valid(char c); - //- Solve the matrix using Gaussian elimination with pivoting - // and return the solution - template - void solve(Field& psi, const Field& source) const; + //- Is the type a wildcard? + inline bool isWildCard() const; - //- LU decompose the matrix with pivoting - static void LUDecompose - ( - Matrix& matrix, - labelList& pivotIndices - ); + // Member operators - //- LU back-substitution with given source, returning the solution - // in the source - template - static void LUBacksubstitute - ( - const Matrix& luMmatrix, - const labelList& pivotIndices, - Field& source - ); + // Assignment - //- Solve the matrix using LU decomposition with pivoting - // returning the LU form of the matrix and the solution in the source - template - static void LUsolve(Matrix& matrix, Field& source); + inline void operator=(const keyType& s); + + //- Assign from regular expression. + inline void operator=(const string& s); + + inline void operator=(const word& s); + + inline void operator=(const char*); + + + // IOstream operators + + friend Istream& operator>>(Istream& is, keyType& w); + + friend Ostream& operator<<(Ostream& os, const keyType& w); }; @@ -114,9 +133,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#ifdef NoRepository -# include "scalarMatrixTemplates.C" -#endif +#include "keyTypeI.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/strings/keyType/keyTypeI.H b/src/OpenFOAM/primitives/strings/keyType/keyTypeI.H new file mode 100644 index 0000000000..f3785ebbff --- /dev/null +++ b/src/OpenFOAM/primitives/strings/keyType/keyTypeI.H @@ -0,0 +1,132 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 + +\*---------------------------------------------------------------------------*/ + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +//- Construct null +inline Foam::keyType::keyType() +: + word(), + isWildCard_(false) +{} + + +//- Construct as copy +inline Foam::keyType::keyType(const keyType& s) +: + word(s, false), + isWildCard_(s.isWildCard()) +{} + + +//- Construct as copy of word +inline Foam::keyType::keyType(const word& s) +: + word(s, false), + isWildCard_(false) +{} + + +//- Construct as copy of string. Expect it to be regular expression +inline Foam::keyType::keyType(const string& s) +: + word(s, false), + isWildCard_(true) +{} + + +//- Construct as copy of character array +inline Foam::keyType::keyType(const char* s) +: + word(s, false), + isWildCard_(false) +{} + + +//- Construct as copy of std::string +inline Foam::keyType::keyType +( + const std::string& s, + const bool isWildCard +) +: + word(s, false), + isWildCard_(isWildCard) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline bool Foam::keyType::valid(char c) +{ + return c != '"'; +} + + +bool Foam::keyType::isWildCard() const +{ + return isWildCard_; +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline void Foam::keyType::operator=(const keyType& s) +{ + // Bypass checking + string::operator=(s); + isWildCard_ = s.isWildCard(); +} + + +inline void Foam::keyType::operator=(const word& s) +{ + word::operator=(s); + isWildCard_ = false; +} + + +inline void Foam::keyType::operator=(const string& s) +{ + // Bypass checking + string::operator=(s); + isWildCard_ = true; +} + + +inline void Foam::keyType::operator=(const char* s) +{ + // Bypass checking + string::operator=(s); + isWildCard_ = false; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/strings/keyType/keyTypeIO.C b/src/OpenFOAM/primitives/strings/keyType/keyTypeIO.C new file mode 100644 index 0000000000..11232282f9 --- /dev/null +++ b/src/OpenFOAM/primitives/strings/keyType/keyTypeIO.C @@ -0,0 +1,88 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 + +Description + Istream constructor and IOstream operators for word. + +\*---------------------------------------------------------------------------*/ + +#include "keyType.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::keyType::keyType(Istream& is) +: + word() +{ + is >> *this; +} + + +Foam::Istream& Foam::operator>>(Istream& is, keyType& w) +{ + token t(is); + + if (!t.good()) + { + is.setBad(); + return is; + } + + if (t.isWord()) + { + w = t.wordToken(); + } + else if (t.isString()) + { + // Assign from string. Sets regular expression. + w = t.stringToken(); + } + else + { + is.setBad(); + FatalIOErrorIn("operator>>(Istream&, keyType&)", is) + << "wrong token type - expected word or string found " + << t.info() + << exit(FatalIOError); + + return is; + } + + // Check state of IOstream + is.check("Istream& operator>>(Istream&, keyType&)"); + + return is; +} + + +Foam::Ostream& Foam::operator<<(Ostream& os, const keyType& w) +{ + os.write(w); + os.check("Ostream& operator<<(Ostream&, const keyType&)"); + return os; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/strings/word/word.H b/src/OpenFOAM/primitives/strings/word/word.H index 3e61275164..6565c28b4f 100644 --- a/src/OpenFOAM/primitives/strings/word/word.H +++ b/src/OpenFOAM/primitives/strings/word/word.H @@ -87,16 +87,21 @@ public: inline word(const word&); //- Construct as copy of character array - inline word(const char*); + inline word(const char*, const bool doStripInvalid = true); //- Construct as copy with a maximum number of characters - inline word(const char*, const size_type); + inline word + ( + const char*, + const size_type, + const bool doStripInvalid + ); //- Construct as copy of string - inline word(const string&); + inline word(const string&, const bool doStripInvalid = true); //- Construct as copy of std::string - inline word(const std::string&); + inline word(const std::string&, const bool doStripInvalid = true); //- Construct from Istream word(Istream&); diff --git a/src/OpenFOAM/primitives/strings/word/wordI.H b/src/OpenFOAM/primitives/strings/word/wordI.H index c31208fd7b..d71eff2f9a 100644 --- a/src/OpenFOAM/primitives/strings/word/wordI.H +++ b/src/OpenFOAM/primitives/strings/word/wordI.H @@ -65,34 +65,51 @@ inline Foam::word::word() {} -inline Foam::word::word(const string& s) +inline Foam::word::word(const string& s, const bool doStripInvalid) : string(s) { - stripInvalid(); + if (doStripInvalid) + { + stripInvalid(); + } } -inline Foam::word::word(const std::string& s) +inline Foam::word::word(const std::string& s, const bool doStripInvalid) : string(s) { - stripInvalid(); + if (doStripInvalid) + { + stripInvalid(); + } } -inline Foam::word::word(const char* s) +inline Foam::word::word(const char* s, const bool doStripInvalid) : string(s) { - stripInvalid(); + if (doStripInvalid) + { + stripInvalid(); + } } -inline Foam::word::word(const char* s, const size_type n) +inline Foam::word::word +( + const char* s, + const size_type n, + const bool doStripInvalid +) : string(s, n) { - stripInvalid(); + if (doStripInvalid) + { + stripInvalid(); + } } diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C index 9629b53bdc..5f5d89c29a 100644 --- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C @@ -826,6 +826,19 @@ void Foam::refinementSurfaces::findNearestIntersection } } } + + + // Make sure that if hit1 has hit something, hit2 will have at least the + // same point (due to tolerances it might miss its end point) + forAll(hit1, pointI) + { + if (hit1[pointI].hit() && !hit2[pointI].hit()) + { + hit2[pointI] = hit1[pointI]; + surface2[pointI] = surface1[pointI]; + region2[pointI] = region1[pointI]; + } + } } diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 0c3fae1887..25137c58c7 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -38,6 +38,9 @@ fvMeshMapper = fvMesh/fvMeshMapper $(fvMeshMapper)/fvPatchMapper.C $(fvMeshMapper)/fvSurfaceMapper.C +extendedStencil = fvMesh/extendedStencil +$(extendedStencil)/extendedStencil.C + fvPatchFields = fields/fvPatchFields $(fvPatchFields)/fvPatchField/fvPatchFields.C @@ -165,6 +168,8 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C +$(schemes)/quadraticFit/quadraticFit.C +$(schemes)/quadraticFit/quadraticFitData.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C new file mode 100644 index 0000000000..bb4c6c1534 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C @@ -0,0 +1,525 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 "extendedStencil.H" +#include "globalIndex.H" +#include "syncTools.H" +#include "SortableList.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Calculates per face a list of global cell/face indices. +void Foam::extendedStencil::calcFaceStencils +( + const polyMesh& mesh, + const globalIndex& globalNumbering +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + const labelList& own = mesh.faceOwner(); + const labelList& nei = mesh.faceNeighbour(); + + + // Determine neighbouring global cell or boundary face + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelList neiGlobal(nBnd); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + // For coupled faces get the cell on the other side + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]); + faceI++; + } + } + else if (isA(pp)) + { + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = -1; + faceI++; + } + } + else + { + // For noncoupled faces get the boundary face. + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = + globalNumbering.toGlobal(mesh.nCells()+bFaceI); + faceI++; + } + } + } + syncTools::swapBoundaryFaceList(mesh, neiGlobal, false); + + + // Determine cellCells in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList globalCellCells(mesh.nCells()); + forAll(globalCellCells, cellI) + { + const cell& cFaces = mesh.cells()[cellI]; + + labelList& cCells = globalCellCells[cellI]; + + cCells.setSize(cFaces.size()); + + // Collect neighbouring cells/faces + label nNbr = 0; + forAll(cFaces, i) + { + label faceI = cFaces[i]; + + if (mesh.isInternalFace(faceI)) + { + label nbrCellI = own[faceI]; + if (nbrCellI == cellI) + { + nbrCellI = nei[faceI]; + } + cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI); + } + else + { + label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()]; + if (nbrCellI != -1) + { + cCells[nNbr++] = nbrCellI; + } + } + } + cCells.setSize(nNbr); + } + + + // Determine neighbouring global cell Cells + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList neiGlobalCellCells(nBnd); + for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + { + neiGlobalCellCells[faceI-mesh.nInternalFaces()] = + globalCellCells[own[faceI]]; + } + syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false); + + + + // Construct stencil in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + stencil_.setSize(mesh.nFaces()); + + labelHashSet faceStencil; + + for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + label globalNei = globalNumbering.toGlobal(nei[faceI]); + faceStencil.insert(globalNei); + const labelList& neiCCells = globalCellCells[nei[faceI]]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + } + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + // Get the coupled cell + label globalNei = neiGlobal[faceI-mesh.nInternalFaces()]; + faceStencil.insert(globalNei); + // And the neighbours of the coupled cell + const labelList& neiCCells = + neiGlobalCellCells[faceI-mesh.nInternalFaces()]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "coupledface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + else if (!isA(pp)) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "boundaryface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + } +} + + +// Calculates extended stencil. This is per face +// - owner +// - cellCells of owner +// - neighbour +// - cellCells of neighbour +// It comes in two parts: +// - a map which collects/distributes all necessary data in a compact array +// - the stencil (a labelList per face) which is a set of indices into this +// compact array. +// The compact array is laid out as follows: +// - first data for current processor (Pstream::myProcNo()) +// - all cells +// - all boundary faces +// - then per processor +// - all used cells and boundary faces +void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) +{ + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + + // Global numbering for cells and boundary faces + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + globalIndex globalNumbering(mesh.nCells()+nBnd); + + + // Calculate stencil in global cell indices + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + calcFaceStencils(mesh, globalNumbering); + + + // Convert stencil to schedule + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // We now know what information we need from other processors. This needs + // to be converted into what information I need to send as well + // (mapDistribute) + + + // 1. Construct per processor compact addressing of the global cells + // needed. The ones from the local processor are not included since + // these are always all needed. + List > globalToProc(Pstream::nProcs()); + { + const labelList& procPatchMap = mesh.globalData().procPatchMap(); + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Presize with (as estimate) size of patch to neighbour. + forAll(procPatchMap, procI) + { + if (procPatchMap[procI] != -1) + { + globalToProc[procI].resize + ( + patches[procPatchMap[procI]].size() + ); + } + } + + // Collect all (non-local) globalcells/faces needed. + forAll(stencil_, faceI) + { + const labelList& stencilCells = stencil_[faceI]; + + forAll(stencilCells, i) + { + label globalCellI = stencilCells[i]; + label procI = globalNumbering.whichProcID(stencilCells[i]); + + if (procI != Pstream::myProcNo()) + { + label nCompact = globalToProc[procI].size(); + globalToProc[procI].insert(globalCellI, nCompact); + } + } + } + // Sort global cells needed (not really necessary) + forAll(globalToProc, procI) + { + if (procI != Pstream::myProcNo()) + { + Map