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/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/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/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/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C index ca7e2d3951..ca9905c037 100644 --- a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C +++ b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C @@ -31,6 +31,7 @@ License #include "volFields.H" #include "globalIndex.H" #include "fvMesh.H" +#include "DynamicField.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -65,11 +66,16 @@ Foam::fileName Foam::externalCoupledFunctionObject::baseDir() const Foam::fileName Foam::externalCoupledFunctionObject::groupDir ( const fileName& commsDir, - const word& regionName, + const word& regionGroupName, const wordRe& groupName ) { - fileName result(commsDir/regionName/string::validate(groupName)); + fileName result + ( + commsDir + /regionGroupName + /string::validate(groupName) + ); result.clean(); return result; @@ -126,11 +132,11 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const if (log_) Info<< type() << ": removing all read files" << endl; - forAll(regionNames_, regionI) + forAll(regionGroupNames_, regionI) { - const word& regionName = regionNames_[regionI]; - const fvMesh& mesh = time_.lookupObject(regionName); - const labelList& groups = regionToGroups_[regionName]; + const word& compName = regionGroupNames_[regionI]; + + const labelList& groups = regionToGroups_[compName]; forAll(groups, i) { label groupI = groups[i]; @@ -141,7 +147,7 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const const word& fieldName = groupReadFields_[groupI][fieldI]; rm ( - groupDir(commsDir_, mesh.dbDir(), groupName) + groupDir(commsDir_, compName, groupName) / fieldName + ".in" ); } @@ -159,22 +165,22 @@ void Foam::externalCoupledFunctionObject::removeWriteFiles() const if (log_) Info<< type() << ": removing all write files" << endl; - forAll(regionNames_, regionI) + forAll(regionGroupNames_, regionI) { - const word& regionName = regionNames_[regionI]; - const fvMesh& mesh = time_.lookupObject(regionName); - const labelList& groups = regionToGroups_[regionName]; + const word& compName = regionGroupNames_[regionI]; + + const labelList& groups = regionToGroups_[compName]; forAll(groups, i) { label groupI = groups[i]; const wordRe& groupName = groupNames_[groupI]; - forAll(groupWriteFields_[groupI], fieldI) + forAll(groupReadFields_[groupI], fieldI) { - const word& fieldName = groupWriteFields_[groupI][fieldI]; + const word& fieldName = groupReadFields_[groupI][fieldI]; rm ( - groupDir(commsDir_, mesh.dbDir(), groupName) + groupDir(commsDir_, compName, groupName) / fieldName + ".out" ); } @@ -376,12 +382,21 @@ void Foam::externalCoupledFunctionObject::readLines void Foam::externalCoupledFunctionObject::writeGeometry ( - const fvMesh& mesh, + const UPtrList& meshes, const fileName& commsDir, const wordRe& groupName ) { - fileName dir(groupDir(commsDir, mesh.dbDir(), groupName)); + wordList regionNames(meshes.size()); + forAll(meshes, i) + { + regionNames[i] = meshes[i].dbDir(); + } + + // Make sure meshes are provided in sorted order + checkOrder(regionNames); + + fileName dir(groupDir(commsDir, compositeName(regionNames), groupName)); //if (log_) { @@ -397,99 +412,210 @@ void Foam::externalCoupledFunctionObject::writeGeometry osFacesPtr.reset(new OFstream(dir/"patchFaces")); } - const labelList patchIDs - ( - mesh.boundaryMesh().patchSet - ( - List(1, groupName) - ).sortedToc() - ); - forAll(patchIDs, i) + DynamicList allMeshesFaces; + DynamicField allMeshesPoints; + + forAll(meshes, meshI) { - label patchI = patchIDs[i]; + const fvMesh& mesh = meshes[meshI]; - const polyPatch& p = mesh.boundaryMesh()[patchI]; + const labelList patchIDs + ( + mesh.boundaryMesh().patchSet + ( + List(1, groupName) + ).sortedToc() + ); + + // Count faces + label nFaces = 0; + forAll(patchIDs, i) + { + nFaces += mesh.boundaryMesh()[patchIDs[i]].size(); + } + + // Collect faces + DynamicList