diff --git a/applications/solvers/basic/potentialFoam/createControls.H b/applications/solvers/basic/potentialFoam/createControls.H deleted file mode 100644 index 7015cae02e..0000000000 --- a/applications/solvers/basic/potentialFoam/createControls.H +++ /dev/null @@ -1,9 +0,0 @@ -const dictionary& potentialFlow -( - mesh.solutionDict().subDict("potentialFlow") -); - -const int nNonOrthCorr -( - potentialFlow.lookupOrDefault("nNonOrthogonalCorrectors", 0) -); diff --git a/applications/solvers/basic/potentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/createFields.H index ce912510ec..36c4dcc08a 100644 --- a/applications/solvers/basic/potentialFoam/createFields.H +++ b/applications/solvers/basic/potentialFoam/createFields.H @@ -12,6 +12,7 @@ volVectorField U mesh ); +// Initialise the velocity internal field to zero U = dimensionedVector("0", U.dimensions(), vector::zero); surfaceScalarField phi @@ -34,7 +35,9 @@ if (args.optionFound("initialiseUBCs")) } -// Default name for the pressure field +// Construct a pressure field +// If it is available read it otherwise construct from the velocity BCs +// converting fixed-value BCs to zero-gradient and vice versa. word pName("p"); // Update name of the pressure field from the command-line option @@ -55,6 +58,9 @@ forAll(U.boundaryField(), patchi) } } +// Note that registerObject is false for the pressure field. The pressure +// field in this solver doesn't have a physical value during the solution. +// It shouldn't be looked up and used by sub models or boundary conditions. Info<< "Constructing pressure field " << pName << nl << endl; volScalarField p ( @@ -64,13 +70,28 @@ volScalarField p runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, - IOobject::NO_WRITE + IOobject::NO_WRITE, + false ), mesh, dimensionedScalar(pName, sqr(dimVelocity), 0), pBCTypes ); +label pRefCell = 0; +scalar pRefValue = 0.0; +if (args.optionFound("writep")) +{ + setRefCell + ( + p, + potentialFlow.dict(), + pRefCell, + pRefValue + ); +} + + Info<< "Constructing velocity potential field Phi\n" << endl; volScalarField Phi ( diff --git a/applications/solvers/basic/potentialFoam/potentialFoam.C b/applications/solvers/basic/potentialFoam/potentialFoam.C index 5c56133788..d6dc4b78a5 100644 --- a/applications/solvers/basic/potentialFoam/potentialFoam.C +++ b/applications/solvers/basic/potentialFoam/potentialFoam.C @@ -24,13 +24,64 @@ License Application potentialFoam -Description - Potential flow solver which solves for the velocity potential - from which the flux-field is obtained and velocity field by reconstructing - the flux. +Group + grpBasicSolvers - This application is particularly useful to generate starting fields for - Navier-Stokes codes. +Description + Potential flow solver. + + \heading Solver details + The potential flow solution is typically employed to generate initial fields + for full Navier-Stokes codes. The flow is evolved using the equation: + + \f[ + \laplacian \Phi = \div(\vec{U}) + \f] + + Where: + \vartable + \Phi | Velocity potential [m2/s] + \vec{U} | Velocity [m/s] + \endvartable + + The corresponding pressure field could be calculated from the divergence + of the Euler equation: + + \f[ + \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0 + \f] + + but this generates excessive pressure variation in regions of large + velocity gradient normal to the flow direction. A better option is to + calculate the pressure field corresponding to velocity variation along the + stream-lines: + + \f[ + \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0 + \f] + where the flow direction tensor \f$\vec{F}\f$ is obtained from + \f[ + \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}} + \f] + + \heading Required fields + \plaintable + U | Velocity [m/s] + \endplaintable + + \heading Optional fields + \plaintable + p | Kinematic pressure [m2/s2] + Phi | Velocity potential [m2/s] + | Generated from p (if present) or U if not present + \endplaintable + + \heading Options + \plaintable + -writep | write the Euler pressure + -writePhi | Write the final velocity potential + -initialiseUBCs | Update the velocity boundaries before solving for Phi + \endplaintable \*---------------------------------------------------------------------------*/ @@ -58,19 +109,13 @@ int main(int argc, char *argv[]) argList::addBoolOption ( "writePhi", - "Write the velocity potential field" + "Write the final velocity potential field" ); argList::addBoolOption ( "writep", - "Calculate and write the pressure field" - ); - - argList::addBoolOption - ( - "withFunctionObjects", - "execute functionObjects" + "Calculate and write the Euler pressure field" ); #include "setRootCase.H" @@ -137,7 +182,7 @@ int main(int argc, char *argv[]) Phi.write(); } - // Calculate the pressure field + // Calculate the pressure field from the Euler equation if (args.optionFound("writep")) { Info<< nl << "Calculating approximate pressure field" << endl; diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C index d2f25bb46f..2fac9f25c0 100644 --- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C +++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C @@ -32,6 +32,7 @@ Description #include "argList.H" #include "Time.H" #include "dynamicFvMesh.H" +#include "pimpleControl.H" #include "vtkSurfaceWriter.H" #include "cyclicAMIPolyPatch.H" @@ -129,13 +130,23 @@ int main(int argc, char *argv[]) Info<< "Writing VTK files with weights of AMI patches." << nl << endl; } + pimpleControl pimple(mesh); + + bool moveMeshOuterCorrectors + ( + pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false) + ); + while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << endl; - for (int i = 0; i<2; i++) + while (pimple.loop()) { - mesh.update(); + if (pimple.firstIter() || moveMeshOuterCorrectors) + { + mesh.update(); + } } mesh.checkMesh(true); diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict index 92b1c02e35..88c86bf482 100644 --- a/applications/utilities/postProcessing/sampling/sample/sampleDict +++ b/applications/utilities/postProcessing/sampling/sample/sampleDict @@ -31,18 +31,66 @@ setFormat raw; // dx : DX scalar or vector format // vtk : VTK ascii format // raw : x y z value format for use with e.g. gnuplot 'splot'. +// boundaryData: written in a form that can be used with timeVaryingMapped +// boundary conditions (move to constant/boundaryData) +// starcd : Star-CD format (geometry in .cel/.vrt, data in .usr) +// nastran : Nastran (.nas) format // // Note: -// other formats such as obj, stl, etc can also be written (by proxy) -// but without any values! +// - other formats such as obj, stl, etc can also be written (by proxy) +// but without any values! +// - boundaryData format: typical use: +// - use a surfaces functionObject to sample a patch: +// type surfaces; +// surfaceFormat boundaryData; +// fields ( p ); +// surfaces +// ( +// inlet +// { +// type patch; +// patches (inpletPatch); +// interpolate false; +// } +// ); +// - the boundaryData writer will write postProcessing/surfaces/inlet/ +// - move this to constant/boundaryData/inlet of the destination case +// - use a timeVaryingMappedFixedValue bc to read&interpolate and use +// as fixedValue. +// - boundaryData: 'interpolate false' writes face centres, 'interpolate true' +// writes points. For 2D case the face centres might only be a single column +// so cannot be used in timeVaryingMapped with standard mapping (since +// uses triangulation to do interpolation). + surfaceFormat vtk; -// optionally define extra controls for the output formats +// Optionally define extra controls for the output formats formatOptions { ensight { + // ascii/binary format format ascii; + //collateTimes true; // write single file containing multiple timesteps + // (only for static surfaces) + } + nastran + { + // From OpenFOAM field name to Nastran field name + fields + ( + (U PLOAD4) + (p PLOAD2) + ); + // Optional scale + scale 1.0; + // Optional format + // short/long/free format + // short: scalar field width 8. default. + // long : scalar field width 16 + // free : comma separated, field width according to writePrecision + // setting + format short; // short/long/free } } @@ -143,9 +191,19 @@ sets type patchSeed; axis xyz; patches (".*Wall.*"); - // Number of points to seed. Divided amongst all processors according - // to fraction of patches they hold. + + // Subset patch faces by: + + // 1. Number of points to seed. Divided amongst all processors + // according to fraction of patches they hold. maxPoints 100; + + // 2. Specified set of locations. This selects for every the specified + // point the nearest patch face. (in addition the number of points + // is also truncated by the maxPoints setting) + // The difference with patchCloud is that this selects patch + // face centres, not an arbitrary location on the face. + points ((0.049 0.099 0.005)(0.051 0.054 0.005)); } ); @@ -249,8 +307,10 @@ surfaces isoValue 0.5; interpolate true; - //zone ABC; // Optional: zone only - //exposedPatchName fixedWalls; // Optional: zone only + // zone ABC; // Optional: zone only + // exposedPatchName fixedWalls; // Optional: zone only + + // bounds (1 1 1)(2 2 2); // Optional: limit extent // regularise false; // Optional: do not simplify // mergeTol 1e-10; // Optional: fraction of mesh bounding box @@ -265,6 +325,9 @@ surfaces isoValue 0.5; interpolate false; regularise false; // do not simplify + + // bounds (1 1 1)(2 2 2); // Optional: limit extent + // mergeTol 1e-10; // Optional: fraction of mesh bounding box // to merge points (default=1e-6) } @@ -281,8 +344,10 @@ surfaces } interpolate true; - //zone ABC; // Optional: zone only - //exposedPatchName fixedWalls; // Optional: zone only + // zone ABC; // Optional: zone only + // exposedPatchName fixedWalls; // Optional: zone only + + // bounds (1 1 1)(2 2 2); // Optional: limit extent // regularise false; // Optional: do not simplify // mergeTol 1e-10; // Optional: fraction of mesh bounding box @@ -305,6 +370,8 @@ surfaces // of isoSurfaceCell interpolate false; regularise false; // Optional: do not simplify + // bounds (1 1 1)(2 2 2); // Optional: limit extent + // mergeTol 1e-10; // Optional: fraction of mesh bounding box // to merge points (default=1e-6) } diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H new file mode 100644 index 0000000000..9f57e35a40 --- /dev/null +++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H @@ -0,0 +1,10 @@ +// +// addRegionOption.H +// ~~~~~~~~~~~~~~~~~ + + Foam::argList::addOption + ( + "regions", + "(name1 .. nameN)", + "specify alternative mesh regions" + ); diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C index 51818b3787..6a8b78c683 100644 --- a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C +++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C @@ -38,6 +38,11 @@ Usage \param -region \ \n Specify an alternative mesh region. + \param -regions (\ \ .. \) \n + Specify alternative mesh regions. The region names will be sorted + alphabetically and a single composite name will be created + \_\.._\ + On execution, the combined patch geometry (points and faces) are output to the communications directory. @@ -59,6 +64,7 @@ SeeAlso int main(int argc, char *argv[]) { #include "addRegionOption.H" + #include "addRegionsOption.H" argList::validArgs.append("patchGroup"); argList::addOption ( @@ -68,16 +74,52 @@ int main(int argc, char *argv[]) ); #include "setRootCase.H" #include "createTime.H" - #include "createNamedMesh.H" - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + wordList regionNames(1, fvMesh::defaultRegion); + if (!args.optionReadIfPresent("region", regionNames[0])) + { + args.optionReadIfPresent("regions", regionNames); + } - const wordRe patchGroup(args[1]); + const wordRe patchGroup(args.argRead(1)); fileName commsDir(runTime.path()/"comms"); args.optionReadIfPresent("commsDir", commsDir); - externalCoupledFunctionObject::writeGeometry(mesh, commsDir, patchGroup); + + // Make sure region names are in canonical order + stableSort(regionNames); + + + PtrList meshes(regionNames.size()); + forAll(regionNames, i) + { + Info<< "Create mesh " << regionNames[i] << " for time = " + << runTime.timeName() << nl << endl; + + meshes.set + ( + i, + new fvMesh + ( + Foam::IOobject + ( + regionNames[i], + runTime.timeName(), + runTime, + Foam::IOobject::MUST_READ + ) + ) + ); + } + + + externalCoupledFunctionObject::writeGeometry + ( + UPtrList(meshes), + commsDir, + patchGroup + ); Info<< "\nEnd\n" << endl; diff --git a/etc/caseDicts/setConstraintTypes b/etc/caseDicts/setConstraintTypes index 6b075c7382..95bce1469d 100644 --- a/etc/caseDicts/setConstraintTypes +++ b/etc/caseDicts/setConstraintTypes @@ -39,11 +39,13 @@ nonuniformTransformCyclic processor { type processor; + value $internalField; } processorCyclic { type processorCyclic; + value $internalField; } symmetryPlane diff --git a/src/OpenFOAM/db/error/IOerror.C b/src/OpenFOAM/db/error/IOerror.C index 08f490c0bd..ee138d80c9 100644 --- a/src/OpenFOAM/db/error/IOerror.C +++ b/src/OpenFOAM/db/error/IOerror.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -174,7 +174,7 @@ void Foam::IOerror::exit(const int) jobInfo.exit(); } - if (abort_) + if (env("FOAM_ABORT")) { abort(); } @@ -215,7 +215,7 @@ void Foam::IOerror::abort() jobInfo.abort(); } - if (abort_) + if (env("FOAM_ABORT")) { Perr<< endl << *this << endl << "\nFOAM aborting (FOAM_ABORT set)\n" << endl; diff --git a/src/OpenFOAM/db/error/error.C b/src/OpenFOAM/db/error/error.C index 91ec556324..f29eafb486 100644 --- a/src/OpenFOAM/db/error/error.C +++ b/src/OpenFOAM/db/error/error.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -40,7 +40,6 @@ Foam::error::error(const string& title) functionName_("unknown"), sourceFileName_("unknown"), sourceFileLineNumber_(0), - abort_(env("FOAM_ABORT")), throwExceptions_(false), messageStreamPtr_(new OStringStream()) { @@ -61,7 +60,6 @@ Foam::error::error(const dictionary& errDict) functionName_(errDict.lookup("functionName")), sourceFileName_(errDict.lookup("sourceFileName")), sourceFileLineNumber_(readLabel(errDict.lookup("sourceFileLineNumber"))), - abort_(env("FOAM_ABORT")), throwExceptions_(false), messageStreamPtr_(new OStringStream()) { @@ -83,7 +81,6 @@ Foam::error::error(const error& err) functionName_(err.functionName_), sourceFileName_(err.sourceFileName_), sourceFileLineNumber_(err.sourceFileLineNumber_), - abort_(err.abort_), throwExceptions_(err.throwExceptions_), messageStreamPtr_(new OStringStream(*err.messageStreamPtr_)) { @@ -173,7 +170,7 @@ void Foam::error::exit(const int errNo) jobInfo.exit(); } - if (abort_) + if (env("FOAM_ABORT")) { abort(); } @@ -214,7 +211,7 @@ void Foam::error::abort() jobInfo.abort(); } - if (abort_) + if (env("FOAM_ABORT")) { Perr<< endl << *this << endl << "\nFOAM aborting (FOAM_ABORT set)\n" << endl; diff --git a/src/OpenFOAM/db/error/error.H b/src/OpenFOAM/db/error/error.H index 97a08b175c..5626ba4db0 100644 --- a/src/OpenFOAM/db/error/error.H +++ b/src/OpenFOAM/db/error/error.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -78,11 +78,10 @@ protected: string sourceFileName_; label sourceFileLineNumber_; - bool abort_; - bool throwExceptions_; OStringStream* messageStreamPtr_; + public: // Constructors diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C new file mode 100644 index 0000000000..8085ccd8f0 --- /dev/null +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C @@ -0,0 +1,242 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "triangle.H" +#include "triPoints.H" +#include "plane.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +template +inline void Foam::triangle::triSliceWithPlane +( + const plane& pl, + const triPoints& tri, + AboveOp& aboveOp, + BelowOp& belowOp +) +{ + // distance to cutting plane + FixedList d; + + // determine how many of the points are above the cutting plane + label nPos = 0; + label posI = -1; + label negI = -1; + forAll(tri, i) + { + d[i] = ((tri[i] - pl.refPoint()) & pl.normal()); + + if (d[i] > 0) + { + nPos++; + posI = i; + } + else + { + negI = i; + } + } + + if (nPos == 3) + { + aboveOp(tri); + } + else if (nPos == 2) + { + // point under the plane + label i0 = negI; + + // indices of remaining points + label i1 = d.fcIndex(i0); + label i2 = d.fcIndex(i1); + + // determine the two intersection points + point p01 = planeIntersection(d, tri, i0, i1); + point p02 = planeIntersection(d, tri, i0, i2); + + aboveOp(triPoints(tri[i1], tri[i2], p02)); + aboveOp(triPoints(tri[i1], p02, p01)); + belowOp(triPoints(tri[i0], p01, p02)); + } + else if (nPos == 1) + { + // point above the plane + label i0 = posI; + + // indices of remaining points + label i1 = d.fcIndex(i0); + label i2 = d.fcIndex(i1); + + // determine the two intersection points + point p01 = planeIntersection(d, tri, i0, i1); + point p02 = planeIntersection(d, tri, i0, i2); + + belowOp(triPoints(tri[i1], tri[i2], p02)); + belowOp(triPoints(tri[i1], p02, p01)); + aboveOp(triPoints(tri[i0], p01, p02)); + } + else + { + // All below + belowOp(tri); + } +} + + +template +template +inline void Foam::triangle::sliceWithPlane +( + const plane& pl, + AboveOp& aboveOp, + BelowOp& belowOp +) const +{ + triSliceWithPlane(pl, triPoints(a_, b_, c_), aboveOp, belowOp); +} + + +template +template +inline void Foam::triangle::triangleOverlap +( + const vector& n, + const triangle& tgt, + InsideOp& insideOp, + OutsideOp& outsideOp +) const +{ + // There are two possibilities with this algorithm - we either cut + // the outside triangles with all the edges or not (and keep them + // as disconnected triangles). In the first case + // we cannot do any evaluation short cut so we've chosen not to re-cut + // the outside triangles. + + + triIntersectionList insideTrisA; + label nInsideA = 0; + storeOp insideOpA(insideTrisA, nInsideA); + + triIntersectionList outsideTrisA; + label nOutsideA = 0; + storeOp outsideOpA(outsideTrisA, nOutsideA); + + + const triPoints thisTri(a_, b_, c_); + + + // Cut original triangle with tgt edge 0. + // From *this to insideTrisA, outsideTrisA. + { + scalar s = Foam::mag(tgt.b() - tgt.a()); + const plane pl0(tgt.a(), tgt.b(), tgt.b() + s*n); + triSliceWithPlane(pl0, thisTri, insideOpA, outsideOpA); + } + + // Shortcut if nothing cut + + if (insideOpA.nTris_ == 0) + { + outsideOp(thisTri); + return; + } + + if (outsideOpA.nTris_ == 0) + { + insideOp(thisTri); + return; + } + + + // Cut all triangles with edge 1. + // From insideTrisA to insideTrisB, outsideTrisA + + triIntersectionList insideTrisB; + label nInsideB = 0; + storeOp insideOpB(insideTrisB, nInsideB); + + //triIntersectionList outsideTrisB; + //label nOutsideB = 0; + //storeOp outsideOpB(outsideTrisB, nOutsideB); + + { + scalar s = Foam::mag(tgt.c() - tgt.b()); + const plane pl0(tgt.b(), tgt.c(), tgt.c() + s*n); + + for (label i = 0; i < insideOpA.nTris_; i++) + { + const triPoints& tri = insideOpA.tris_[i]; + triSliceWithPlane(pl0, tri, insideOpB, outsideOpA); + } + + //// Recut outside triangles (not necessary if only interested in + //// intersection properties) + //for (label i = 0; i < outsideOpA.nTris_; i++) + //{ + // const triPoints& tri = outsideOpA.tris_[i]; + // triSliceWithPlane(pl0, tri, outsideOpB, outsideOpB); + //} + } + + + // Cut all triangles with edge 2. + // From insideTrisB to insideTrisA, outsideTrisA + { + scalar s = Foam::mag(tgt.a() - tgt.c()); + const plane pl0(tgt.c(), tgt.a(), tgt.a() + s*n); + + insideOpA.nTris_ = 0; + //outsideOpA.nTris_ = 0; + for (label i = 0; i < insideOpB.nTris_; i++) + { + const triPoints& tri = insideOpB.tris_[i]; + triSliceWithPlane(pl0, tri, insideOpA, outsideOpA); + } + + //// Recut outside triangles (not necessary if only interested in + //// intersection properties) + //for (label i = 0; i < outsideOpB.nTris_; i++) + //{ + // const triPoints& tri = outsideOpB.tris_[i]; + // triSliceWithPlane(pl0, tri, outsideOpA, outsideOpA); + //} + } + + // Transfer from A to argument + for (label i = 0; i < insideOpA.nTris_; i++) + { + insideOp(insideOpA.tris_[i]); + } + + for (label i = 0; i < outsideOpA.nTris_; i++) + { + outsideOp(outsideOpA.tris_[i]); + } +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 6e8030c797..de477d20e0 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -52,6 +52,8 @@ namespace Foam class Istream; class Ostream; +class triPoints; +class plane; // Forward declaration of friend functions and operators @@ -71,6 +73,8 @@ inline Ostream& operator<< const triangle& ); +typedef triangle triPointRef; + /*---------------------------------------------------------------------------*\ Class triangle Declaration @@ -79,27 +83,98 @@ inline Ostream& operator<< template class triangle { +public: + + // Public typedefs + + //- Storage type for triangles originating from intersecting triangle + // with another triangle + typedef FixedList triIntersectionList; + + //- Return types for classify + enum proxType + { + NONE, + POINT, // Close to point + EDGE // Close to edge + }; + + + // Public classes + + // Classes for use in sliceWithPlane. What to do with decomposition + // of triangle. + + //- Dummy + class dummyOp + { + public: + inline void operator()(const triPoints&); + }; + + //- Sum resulting areas + class sumAreaOp + { + public: + scalar area_; + + inline sumAreaOp(); + + inline void operator()(const triPoints&); + }; + + //- Store resulting tris + class storeOp + { + public: + triIntersectionList& tris_; + label& nTris_; + + inline storeOp(triIntersectionList&, label&); + + inline void operator()(const triPoints&); + }; + + +private: + // Private data PointRef a_, b_, c_; + // Private Member Functions + + //- Helper: calculate intersection point + inline static point planeIntersection + ( + const FixedList& d, + const triPoints& t, + const label negI, + const label posI + ); + + //- Helper: slice triangle with plane + template + inline static void triSliceWithPlane + ( + const plane& pl, + const triPoints& tri, + AboveOp& aboveOp, + BelowOp& belowOp + ); + + public: - //- Return types for classify - enum proxType - { - NONE, - POINT, // Close to point - EDGE // Close to edge - }; - - // Constructors //- Construct from three points inline triangle(const Point& a, const Point& b, const Point& c); + //- Construct from three points + inline triangle(const FixedList&); + //- Construct from three points in the list of points // The indices could be from triFace etc. inline triangle @@ -237,6 +312,27 @@ public: pointHit& edgePoint ) const; + //- Decompose triangle into triangles above and below plane + template + inline void sliceWithPlane + ( + const plane& pl, + AboveOp& aboveOp, + BelowOp& belowOp + ) const; + + //- Decompose triangle into triangles inside and outside + // (with respect to user provided normal) other + // triangle. + template + inline void triangleOverlap + ( + const vector& n, + const triangle& tri, + InsideOp& insideOp, + OutsideOp& outsideOp + ) const; + // IOstream operators @@ -264,6 +360,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "triangle.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index fbc5a14286..e3f8993e2e 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,6 +25,7 @@ License #include "IOstreams.H" #include "pointHit.H" +#include "triPoints.H" #include "mathematicalConstants.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -43,6 +44,18 @@ inline Foam::triangle::triangle {} +template +inline Foam::triangle::triangle +( + const FixedList& tri +) +: + a_(tri[0]), + b_(tri[1]), + c_(tri[2]) +{} + + template inline Foam::triangle::triangle ( @@ -804,6 +817,66 @@ inline Foam::pointHit Foam::triangle::nearestPoint } +template +inline void Foam::triangle::dummyOp::operator() +( + const triPoints& +) +{} + + +template +inline Foam::triangle::sumAreaOp::sumAreaOp() +: + area_(0.0) +{} + + +template +inline void Foam::triangle::sumAreaOp::operator() +( + const triPoints& tri +) +{ + area_ += triangle(tri).mag(); +} + + +template +inline Foam::triangle::storeOp::storeOp +( + triIntersectionList& tris, + label& nTris +) +: + tris_(tris), + nTris_(nTris) +{} + + +template +inline void Foam::triangle::storeOp::operator() +( + const triPoints& tri +) +{ + tris_[nTris_++] = tri; +} + + +template +inline Foam::point Foam::triangle::planeIntersection +( + const FixedList& d, + const triPoints& t, + const label negI, + const label posI +) +{ + return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]); +} + + // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // template diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C index d546ff86bb..dbef33a6ad 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C @@ -514,12 +514,31 @@ void Foam::motionSmootherAlgo::setDisplacement const labelList& cppMeshPoints = mesh.globalData().coupledPatch().meshPoints(); - forAll(cppMeshPoints, i) + const labelList& ppMeshPoints = pp.meshPoints(); + + // Knock out displacement on points which are not on pp but are coupled + // to them since we want 'proper' values from displacement to take + // precedence. { - displacement[cppMeshPoints[i]] = vector::zero; + PackedBoolList isPatchPoint(mesh.nPoints()); + isPatchPoint.set(ppMeshPoints); + syncTools::syncPointList + ( + mesh, + isPatchPoint, + maxEqOp(), + 0 + ); + forAll(cppMeshPoints, i) + { + label pointI = cppMeshPoints[i]; + if (isPatchPoint[pointI]) + { + displacement[pointI] = vector::zero; + } + } } - const labelList& ppMeshPoints = pp.meshPoints(); // Set internal point data from displacement on combined patch points. forAll(ppMeshPoints, patchPointI) diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H index 724658012e..4cc17efdcd 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H @@ -56,7 +56,8 @@ Note and/or edges but no faces of pp). Hence we have to be careful when e.g. synchronising displacements that the value from the processor which has faces of pp get priority. This is currently handled in setDisplacement - by resetting the internal displacement to zero before doing anything + by resetting the internal displacement to zero on coupled points + that are coupled to patch points before doing anything else. The combine operator used will give preference to non-zero values. diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C index a4148b950b..b4117c424a 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -48,12 +48,9 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio return 1.0; } - //scalar totalOverFirst = totalThickness/firstLayerThickess; - - const label maxIters = 10; + const label maxIters = 20; const scalar tol = 1e-8; - if (mag(n-totalOverFirst) < tol) { return 1.0; @@ -74,8 +71,6 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio maxR = totalOverFirst/(n - 1); } - //Info<< "Solution bounds = (" << minR << ", " << maxR << ")" << nl << endl; - // Starting guess scalar r = 0.5*(minR + maxR); @@ -85,14 +80,9 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio const scalar fx = pow(r, n) - totalOverFirst*r - (1 - totalOverFirst); const scalar dfx = n*pow(r, n - 1) - totalOverFirst; - r -= fx/dfx; - const scalar error = mag(r - prevr); - - //Info<< i << " " << r << " Error = " << error << endl; - - if (error < tol) + if (mag(r - prevr) < tol) { break; } @@ -103,7 +93,6 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -// Construct from dictionary Foam::layerParameters::layerParameters ( const dictionary& dict, @@ -417,9 +406,9 @@ Foam::scalar Foam::layerParameters::layerThickness } else { - return firstLayerThickess * - (1.0 - pow(expansionRatio, nLayers)) - / (1.0 - expansionRatio); + return firstLayerThickess + *(1.0 - pow(expansionRatio, nLayers)) + /(1.0 - expansionRatio); } } break; @@ -433,9 +422,9 @@ Foam::scalar Foam::layerParameters::layerThickness else { scalar invExpansion = 1.0 / expansionRatio; - return finalLayerThickess * - (1.0 - pow(invExpansion, nLayers)) - / (1.0 - invExpansion); + return finalLayerThickess + *(1.0 - pow(invExpansion, nLayers)) + /(1.0 - invExpansion); } } break; @@ -484,7 +473,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio { return 1.0 - / layerExpansionRatio + /layerExpansionRatio ( nLayers, totalThickness/finalLayerThickess @@ -577,8 +566,8 @@ Foam::scalar Foam::layerParameters::finalLayerThicknessRatio { return pow(expansionRatio, nLayers - 1) - * (1.0 - expansionRatio) - / (1.0 - pow(expansionRatio, nLayers)); + *(1.0 - expansionRatio) + /(1.0 - pow(expansionRatio, nLayers)); } } else diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H index 505fb8b07b..f617916f47 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H @@ -2,8 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -333,7 +333,6 @@ public: const label nLayers, const scalar expansionRatio ) const; - }; diff --git a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H index 9fea06cf2c..b832f7e2be 100644 --- a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H +++ b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H @@ -33,8 +33,8 @@ SeeAlso \*---------------------------------------------------------------------------*/ -#ifndef plane_H -#define plane_H +#ifndef planeExtrusion_H +#define planeExtrusion_H #include "linearNormal.H" diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C index 7a862ee9da..84028ebef5 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C @@ -44,7 +44,7 @@ const Foam::scalar Foam::cyclicACMIPolyPatch::tolerance_ = 1e-6; void Foam::cyclicACMIPolyPatch::initPatchFaceAreas() const { - if (!empty() && faceAreas0_.empty()) + if (!empty() && (faceAreas0_.empty() || boundaryMesh().mesh().moving())) { faceAreas0_ = faceAreas(); } @@ -52,9 +52,13 @@ void Foam::cyclicACMIPolyPatch::initPatchFaceAreas() const const cyclicACMIPolyPatch& nbrACMI = refCast(this->neighbPatch()); - if (!nbrACMI.empty() && nbrACMI.faceAreas0().empty()) + if + ( + !nbrACMI.empty() + && (nbrACMI.faceAreas0().empty() || boundaryMesh().mesh().moving()) + ) { - nbrACMI.initPatchFaceAreas(); + nbrACMI.faceAreas0_ = nbrACMI.faceAreas(); } } @@ -136,11 +140,13 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const void Foam::cyclicACMIPolyPatch::initGeometry(PstreamBuffers& pBufs) { - // Initialise the AMI so that base geometry (e.g. cell volumes) are - // correctly evaluated - resetAMI(); - + // Note: cyclicAMIPolyPatch clears AMI so do first cyclicAMIPolyPatch::initGeometry(pBufs); + + // Initialise the AMI so that base geometry (e.g. cell volumes) are + // correctly evaluated before e.g. any of the processor patches gets + // hit (since uses cell volumes in its initGeometry) + resetAMI(); } @@ -156,7 +162,13 @@ void Foam::cyclicACMIPolyPatch::initMovePoints const pointField& p ) { + // Note: cyclicAMIPolyPatch clears AMI so do first cyclicAMIPolyPatch::initMovePoints(pBufs, p); + + // Initialise the AMI so that base geometry (e.g. cell volumes) are + // correctly evaluated before e.g. any of the processor patches gets + // hit (since uses cell volumes in its initGeometry) + resetAMI(); } diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C index 231a2e0ba7..a3a3707e01 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,7 +37,7 @@ Foam::tmp > Foam::cyclicACMIPolyPatch::interpolate if (owner()) { - const scalarField& w = srcMask_; + const scalarField& w = AMI().srcWeightsSum(); tmp > interpField(AMI().interpolateToSource(fldCouple)); @@ -45,7 +45,7 @@ Foam::tmp > Foam::cyclicACMIPolyPatch::interpolate } else { - const scalarField& w = neighbPatch().tgtMask(); + const scalarField& w = neighbPatch().AMI().tgtWeightsSum(); tmp > interpField ( @@ -82,16 +82,18 @@ void Foam::cyclicACMIPolyPatch::interpolate if (owner()) { - const scalarField& w = srcMask_; + const scalarField& w = AMI().srcWeightsSum(); AMI().interpolateToSource(fldCouple, cop, result); + result = result + (1.0 - w)*fldNonOverlap; } else { - const scalarField& w = neighbPatch().tgtMask(); + const scalarField& w = neighbPatch().AMI().tgtWeightsSum(); neighbPatch().AMI().interpolateToTarget(fldCouple, cop, result); + result = result + (1.0 - w)*fldNonOverlap; } } diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C index f58e99ce0e..b8f9b3b2c0 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C @@ -24,9 +24,13 @@ License \*---------------------------------------------------------------------------*/ #include "cyclicAMIPolyPatch.H" +#include "transformField.H" #include "SubField.H" +#include "polyMesh.H" #include "Time.H" #include "addToRunTimeSelectionTable.H" +#include "faceAreaIntersect.H" +#include "ops.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -172,18 +176,19 @@ void Foam::cyclicAMIPolyPatch::calcTransforms scalar errorPos = mag(transformedAreaPos + area0); scalar errorNeg = mag(transformedAreaNeg + area0); - if (errorPos < errorNeg) - { - revT = revTPos; - } - else + scalar scaledErrorPos = errorPos/(mag(area0) + ROOTVSMALL); + scalar scaledErrorNeg = errorNeg/(mag(area0) + ROOTVSMALL); + + // One of the errors should be (close to) zero. If this is + // the reverse transformation flip the rotation angle. + revT = revTPos; + if (errorPos > errorNeg && scaledErrorNeg <= matchTolerance()) { revT = revTNeg; rotationAngle_ *= -1; } - scalar areaError = - min(errorPos, errorNeg)/(mag(area0) + ROOTVSMALL); + scalar areaError = min(scaledErrorPos, scaledErrorNeg); if (areaError > matchTolerance()) { @@ -406,6 +411,9 @@ void Foam::cyclicAMIPolyPatch::resetAMI void Foam::cyclicAMIPolyPatch::initGeometry(PstreamBuffers& pBufs) { + // The AMI is no longer valid. Leave it up to demand-driven calculation + AMIPtr_.clear(); + polyPatch::initGeometry(pBufs); } @@ -431,6 +439,9 @@ void Foam::cyclicAMIPolyPatch::initMovePoints const pointField& p ) { + // The AMI is no longer valid. Leave it up to demand-driven calculation + AMIPtr_.clear(); + polyPatch::initMovePoints(pBufs, p); // See below. Clear out any local geometry @@ -447,19 +458,15 @@ void Foam::cyclicAMIPolyPatch::movePoints polyPatch::movePoints(pBufs, p); calcTransforms(); - - // Note: resetAMI is called whilst in geometry update. So the slave - // side might not have reached 'movePoints'. Is explicitly handled by - // - clearing geometry of neighbour inside initMovePoints - // - not using localPoints() inside resetAMI - resetAMI(); } void Foam::cyclicAMIPolyPatch::initUpdateMesh(PstreamBuffers& pBufs) { - polyPatch::initUpdateMesh(pBufs); + // The AMI is no longer valid. Leave it up to demand-driven calculation AMIPtr_.clear(); + + polyPatch::initUpdateMesh(pBufs); } diff --git a/src/meshTools/regionSplit/regionSplit.C b/src/meshTools/regionSplit/regionSplit.C index 2583c2ddfe..1a8de0fa15 100644 --- a/src/meshTools/regionSplit/regionSplit.C +++ b/src/meshTools/regionSplit/regionSplit.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,6 +28,8 @@ License #include "processorPolyPatch.H" #include "globalIndex.H" #include "syncTools.H" +#include "FaceCellWave.H" +#include "minData.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -40,340 +42,93 @@ defineTypeNameAndDebug(regionSplit, 0); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -// Handle (non-processor) coupled faces. -void Foam::regionSplit::transferCoupledFaceRegion -( - const label faceI, - const label otherFaceI, - - labelList& faceRegion, - DynamicList