diff --git a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H index 371d95a0632d1a0de12d9b91109eef9de8eb9c4d..653da9c46d48e8f4a6937d1f9251ec56681d37bb 100644 --- a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H +++ b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -92,6 +92,14 @@ public: const bool doCollapse = false ); + //- Construct from components + inline cellShape + ( + const word& model, + const labelList&, + const bool doCollapse = false + ); + //- Construct from Istream inline cellShape(Istream& is); diff --git a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H index 5120cd4dfa51234fc24eee1dc019a24430b90ad2..42055e04c9229190799ceb06d9b0149972ee7e26 100644 --- a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H +++ b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H @@ -25,6 +25,7 @@ License #include "Istream.H" #include "cell.H" +#include "cellModeller.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -51,6 +52,23 @@ inline Foam::cellShape::cellShape } +inline Foam::cellShape::cellShape +( + const word& model, + const labelList& l, + const bool doCollapse +) +: + labelList(l), + m(cellModeller::lookup(model)) +{ + if (doCollapse) + { + collapse(); + } +} + + inline Foam::cellShape::cellShape(Istream& is) { is >> *this; diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index 4d1da72068ac7f76200b87edacdd1c8fced755e4..f552dec8f9d72c462ec96552b552f18698d29192 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -1,6 +1,7 @@ blockVertices/blockVertex/blockVertex.C blockVertices/pointVertex/pointVertex.C blockVertices/projectVertex/projectVertex.C +blockVertices/namedVertex/namedVertex.C blockEdges/blockEdge/blockEdge.C blockEdges/lineDivide/lineDivide.C @@ -13,6 +14,7 @@ blockEdges/BSplineEdge/BSplineEdge.C blockEdges/splineEdge/CatmullRomSpline.C blockEdges/splineEdge/splineEdge.C blockEdges/projectEdge/projectEdge.C +blockEdges/projectCurveEdge/projectCurveEdge.C blockFaces/blockFace/blockFace.C blockFaces/projectFace/projectFace.C @@ -23,8 +25,9 @@ gradingDescriptor/gradingDescriptors.C blockDescriptor/blockDescriptor.C blockDescriptor/blockDescriptorEdges.C -block/block.C -block/blockCreate.C +blocks/block/block.C +blocks/block/blockCreate.C +blocks/namedBlock/namedBlock.C blockMesh/blockMesh.C blockMesh/blockMeshCreate.C @@ -33,4 +36,6 @@ blockMesh/blockMeshCheck.C blockMesh/blockMeshMerge.C blockMesh/blockMeshMergeFast.C +searchableCurve/searchableCurve.C + LIB = $(FOAM_LIBBIN)/libblockMesh diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options index bf6ca71abfb4c715b09e5ed1679f365702c55784..5467dc3ce38635ce1b666c5263855f795f66e302 100644 --- a/src/mesh/blockMesh/Make/options +++ b/src/mesh/blockMesh/Make/options @@ -1,9 +1,11 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/fileFormats/lnInclude \ + -I$(LIB_SRC)/edgeMesh/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude LIB_LIBS = \ -lmeshTools \ -lfileFormats \ + -ledgeMesh \ -lsurfMesh diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C index fb9a413013b96eab4cb401f21335a6b16d896027..ddcb9466e048a2c14256a5e93178ab33e9654b65 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C @@ -128,6 +128,90 @@ void Foam::blockDescriptor::findCurvedFaces() } +void Foam::blockDescriptor::read +( + Istream& is, + label& val, + const dictionary& dict +) +{ + token t(is); + if (t.isLabel()) + { + val = t.labelToken(); + } + else if (t.isWord()) + { + const word& varName = t.wordToken(); + const entry* ePtr = dict.lookupScopedEntryPtr + ( + varName, + true, + true + ); + if (ePtr) + { + // Read as label + val = Foam::readLabel(ePtr->stream()); + } + else + { + FatalIOErrorInFunction(is) + << "Undefined variable " + << varName << ". Valid variables are " << dict + << exit(FatalIOError); + } + } + else + { + FatalIOErrorInFunction(is) + << "Illegal token " << t.info() + << " when trying to read label" + << exit(FatalIOError); + } + + is.fatalCheck + ( + "operator>>(Istream&, List<T>&) : reading entry" + ); +} + + +Foam::label Foam::blockDescriptor::read +( + Istream& is, + const dictionary& dict +) +{ + label val; + read(is, val, dict); + return val; +} + + +void Foam::blockDescriptor::write +( + Ostream& os, + const label val, + const dictionary& dict +) +{ + forAllConstIter(dictionary, dict, iter) + { + if (iter().isStream()) + { + label keyVal(Foam::readLabel(iter().stream())); + if (keyVal == val) + { + os << iter().keyword(); + return; + } + } + } + os << val; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::blockDescriptor::blockDescriptor @@ -164,6 +248,8 @@ Foam::blockDescriptor::blockDescriptor Foam::blockDescriptor::blockDescriptor ( + const dictionary& dict, + const label index, const pointField& vertices, const blockEdgeList& edges, const blockFaceList& faces, @@ -173,13 +259,24 @@ Foam::blockDescriptor::blockDescriptor vertices_(vertices), edges_(edges), faces_(faces), - blockShape_(is), density_(), expand_(12, gradingDescriptors()), zoneName_(), curvedFaces_(-1), nCurvedFaces_(0) { + // Read cell model and list of vertices (potentially with variables) + word model(is); + blockShape_ = cellShape + ( + model, + read<label> + ( + is, + dict.subOrEmptyDict("namedVertices") + ) + ); + // Examine next token token t(is); diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H index b0f7b1412cfef83aa8dca994b71f5df70864ceee..45c5e755614df75730e4e9d7e90584b6a121a59b 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H @@ -152,6 +152,8 @@ public: //- Construct from Istream blockDescriptor ( + const dictionary& dict, + const label index, const pointField& vertices, const blockEdgeList&, const blockFaceList&, @@ -161,6 +163,9 @@ public: // Member Functions + //- Reference to point field defining the block mesh + inline const pointField& vertices() const; + //- Return reference to the list of curved faces inline const blockFaceList& faces() const; @@ -233,6 +238,23 @@ public: // to lie on the faces of the block void correctFacePoints(FixedList<pointField, 6>&) const; + //- In-place read with dictionary lookup + static void read(Istream&, label&, const dictionary&); + + //- In-place read with dictionary lookup + template<class T> + static void read(Istream&, List<T>&, const dictionary&); + + //- Return-read with dictionary lookup + static label read(Istream&, const dictionary&); + + //- Return-read with dictionary lookup + template<class T> + static List<T> read(Istream& is, const dictionary&); + + //- Write with dictionary lookup + static void write(Ostream&, const label, const dictionary&); + // IOstream Operators @@ -248,6 +270,10 @@ public: #include "blockDescriptorI.H" +#ifdef NoRepository + #include "blockDescriptorTemplates.C" +#endif + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H index 823860a820d276c6693f745195d24c17bc17900a..13ec39c2e67f9798c2cbbd763439406c4012e742 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H @@ -25,6 +25,12 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline const Foam::pointField& Foam::blockDescriptor::vertices() const +{ + return vertices_; +} + + inline const Foam::blockFaceList& Foam::blockDescriptor::faces() const { return faces_; diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..aacd20d8bd32bb56dd23e95beea28de75a79bd24 --- /dev/null +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C @@ -0,0 +1,114 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class T> +void Foam::blockDescriptor::read +( + Istream& is, + List<T>& L, + const dictionary& dict +) +{ + token firstToken(is); + + if (firstToken.isLabel()) + { + label s = firstToken.labelToken(); + + // Set list length to that read + L.setSize(s); + + // Read list contents depending on data format + + // Read beginning of contents + char delimiter = is.readBeginList("List"); + + if (s) + { + if (delimiter == token::BEGIN_LIST) + { + for (label i=0; i<s; i++) + { + read(is, L[i], dict); + } + } + } + + // Read end of contents + is.readEndList("List"); + } + else if (firstToken.isPunctuation()) + { + if (firstToken.pToken() != token::BEGIN_LIST) + { + FatalIOErrorInFunction(is) + << "incorrect first token, expected '(', found " + << firstToken.info() + << exit(FatalIOError); + } + + SLList<T> sll; + + while (true) + { + token t(is); + if (t.isPunctuation() && t.pToken() == token::END_LIST) + { + break; + } + is.putBack(t); + T elem; + read(is, elem, dict); + sll.append(elem); + } + + // Convert the singly-linked list to this list + L = sll; + } + else + { + FatalIOErrorInFunction(is) + << "incorrect first token, expected <int> or '(', found " + << firstToken.info() + << exit(FatalIOError); + } +} + + +template<class T> +Foam::List<T> Foam::blockDescriptor::read +( + Istream& is, + const dictionary& dict +) +{ + List<T> L; + read(is, L, dict); + return L; +} + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C index 4e16c18264c06d3925c628910180db9efefa51f1..dc62dbc3e4a47fb8a41b197c5adc7f633113e6d5 100644 --- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C +++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C @@ -62,12 +62,14 @@ Foam::blockEdges::BSplineEdge::BSplineEdge Foam::blockEdges::BSplineEdge::BSplineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ) : - blockEdge(points, is), + blockEdge(dict, index, points, is), BSpline(appendEndPoints(points, start_, end_, pointField(is))) { token t(is); diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H index 81c06d63a79041d853429d2ae34fe04de6b748df..72cbeeab7ac9d31c2c05236c3d2aa35280f31b04 100644 --- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H +++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H @@ -83,6 +83,8 @@ public: //- Construct from Istream, setting pointsList BSplineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField&, Istream& diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C index 5d012e0531c35ebe7adae77608e9e55afc5e0443..793cc27b2f58e3eda49e5245cb502aa3b685a1ca 100644 --- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C +++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C @@ -123,12 +123,14 @@ Foam::blockEdges::arcEdge::arcEdge Foam::blockEdges::arcEdge::arcEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ) : - blockEdge(points, is), + blockEdge(dict, index, points, is), p1_(points_[start_]), p2_(is), p3_(points_[end_]), diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H index b7a71ff6d102fb56a66e9f3239832c3d0bc66be2..9d41e8cceb7f07a2d48124facecfab9b5619b9ef 100644 --- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H +++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H @@ -91,6 +91,8 @@ public: //- Construct from Istream setting pointsList arcEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C index 4ed2647a741782235ee66daf06059a75a9adb30f..0ef5fd0408e5aa2ee5dfd945b050140fd1de2926 100644 --- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C +++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "blockEdge.H" +#include "blockDescriptor.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -51,13 +52,15 @@ Foam::blockEdge::blockEdge Foam::blockEdge::blockEdge ( + const dictionary& dict, + const label index, const pointField& points, Istream& is ) : points_(points), - start_(readLabel(is)), - end_(readLabel(is)) + start_(blockDescriptor::read(is, dict.subOrEmptyDict("namedVertices"))), + end_(blockDescriptor::read(is, dict.subOrEmptyDict("namedVertices"))) {} @@ -70,6 +73,8 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::clone() const Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is @@ -95,7 +100,7 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New << abort(FatalError); } - return autoPtr<blockEdge>(cstrIter()(geometry, points, is)); + return autoPtr<blockEdge>(cstrIter()(dict, index, geometry, points, is)); } @@ -139,6 +144,25 @@ Foam::blockEdge::position(const scalarList& lambdas) const } +void Foam::blockEdge::write(Ostream& os, const dictionary& d) const +{ + const dictionary* varDictPtr = d.subDictPtr("namedVertices"); + if (varDictPtr) + { + const dictionary& varDict = *varDictPtr; + + blockDescriptor::write(os, start_, varDict); + os << tab; + blockDescriptor::write(os, end_, varDict); + os << endl; + } + else + { + os << start_ << tab << end_ << endl; + } +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const blockEdge& p) diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H index e6e3196ea94a3981245aff19e4368af12e1550ce..bde77adbc29b4feee21aa483a361e84e30a27153 100644 --- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H +++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H @@ -92,11 +92,13 @@ public: blockEdge, Istream, ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ), - (geometry, points, is) + (dict, index, geometry, points, is) ); @@ -113,6 +115,8 @@ public: //- Construct from Istream setting pointsList blockEdge ( + const dictionary& dict, + const label index, const pointField&, Istream& ); @@ -123,6 +127,8 @@ public: //- New function which constructs and returns pointer to a blockEdge static autoPtr<blockEdge> New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField&, Istream& @@ -132,20 +138,29 @@ public: // PtrLists of blockEdge class iNew { + const dictionary& dict_; const searchableSurfaces& geometry_; const pointField& points_; + mutable label index_; public: - iNew(const searchableSurfaces& geometry, const pointField& points) + iNew + ( + const dictionary& dict, + const searchableSurfaces& geometry, + const pointField& points + ) : + dict_(dict), geometry_(geometry), - points_(points) + points_(points), + index_(0) {} autoPtr<blockEdge> operator()(Istream& is) const { - return blockEdge::New(geometry_, points_, is); + return blockEdge::New(dict_, index_++, geometry_, points_, is); } }; @@ -195,6 +210,9 @@ public: //- Return the length of the curve virtual scalar length() const = 0; + //- Write edge with variable backsubstitution + void write(Ostream&, const dictionary&) const; + // Ostream operator diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C index 45b4434850a1f3cb5964e86836c08aa143e2e25a..9206151f56223b0e7809cff921f36a4e7aeaca6e 100644 --- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C +++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C @@ -53,12 +53,14 @@ Foam::blockEdges::lineEdge::lineEdge Foam::blockEdges::lineEdge::lineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ) : - blockEdge(points, is) + blockEdge(dict, index, points, is) {} diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H index e70ffd10bf79477159eaa3e3ed15729963361200..c5266b28c3d2c49d4e7c70c9a4a3c8894db25d94 100644 --- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H +++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H @@ -68,6 +68,8 @@ public: //- Construct from Istream with a pointField lineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField&, Istream& diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C index 1396636a66f62ebe1ed0be3d969e712e52a169d3..63e8fb4d00be11e19a900ed9dd7a542fad7b1499 100644 --- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C +++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C @@ -55,12 +55,14 @@ Foam::blockEdges::polyLineEdge::polyLineEdge Foam::blockEdges::polyLineEdge::polyLineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& ps, Istream& is ) : - blockEdge(ps, is), + blockEdge(dict, index, ps, is), polyLine(appendEndPoints(ps, start_, end_, pointField(is))) {} diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H index adcd7fd1b47c641467ad89d1c2513da007bbef9d..62e42ca31c302916bec0d0ebd0708f0141dbf69c 100644 --- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H +++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H @@ -83,6 +83,8 @@ public: //- Construct from Istream polyLineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField&, Istream& diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C new file mode 100644 index 0000000000000000000000000000000000000000..0e0166fd4463c3493b0a56ee583e04f7e749afff --- /dev/null +++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C @@ -0,0 +1,255 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "searchableSurfacesQueries.H" +#include "projectCurveEdge.H" +#include "unitConversion.H" +#include "addToRunTimeSelectionTable.H" +#include "pointConstraint.H" +#include "OBJstream.H" +#include "linearInterpolationWeights.H" +#include "searchableCurve.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(projectCurveEdge, 0); + addToRunTimeSelectionTable(blockEdge, projectCurveEdge, Istream); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::projectCurveEdge::projectCurveEdge +( + const dictionary& dict, + const label index, + const searchableSurfaces& geometry, + const pointField& points, + Istream& is +) +: + blockEdge(dict, index, points, is), + geometry_(geometry) +{ + wordList names(is); + surfaces_.setSize(names.size()); + forAll(names, i) + { + surfaces_[i] = geometry_.findSurfaceID(names[i]); + + if (surfaces_[i] == -1) + { + FatalIOErrorInFunction(is) + << "Cannot find surface " << names[i] << " in geometry" + << exit(FatalIOError); + } + + if (isA<searchableCurve>(geometry_[surfaces_[i]])) + { + Info<< type() << " : Using curved surface " + << geometry_[surfaces_[i]].name() + << " to predict starting points." << endl; + } + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp<Foam::pointField> +Foam::projectCurveEdge::position(const scalarList& lambdas) const +{ + // For debugging to tag the output + static label eIter = 0; + + autoPtr<OBJstream> debugStr; + if (debug) + { + debugStr.reset + ( + new OBJstream("projectCurveEdge_" + Foam::name(eIter++) + ".obj") + ); + Info<< "Writing lines from straight-line start points" + << " to projected points to " << debugStr().name() << endl; + } + + + tmp<pointField> tpoints(new pointField(lambdas.size())); + pointField& points = tpoints.ref(); + + const point& startPt = points_[start_]; + const point& endPt = points_[end_]; + const vector d = endPt-startPt; + + // Initial guess + forAll(lambdas, i) + { + points[i] = startPt+lambdas[i]*d; + } + + // Use special interpolation to keep initial guess on same position on + // surface + forAll(surfaces_, i) + { + if (isA<searchableCurve>(geometry_[surfaces_[i]])) + { + const searchableCurve& s = + refCast<const searchableCurve>(geometry_[surfaces_[i]]); + List<pointIndexHit> nearInfo; + s.findNearest + ( + points[0], + points.last(), + scalarField(lambdas), + scalarField(points.size(), magSqr(d)), + nearInfo + ); + forAll(nearInfo, i) + { + if (nearInfo[i].hit()) + { + points[i] = nearInfo[i].hitPoint(); + } + } + + break; + } + } + + + + // Upper limit for number of iterations + const label maxIter = 10; + // Residual tolerance + const scalar relTol = 0.1; + const scalar absTol = 1e-4; + + scalar initialResidual = 0.0; + + for (label iter = 0; iter < maxIter; iter++) + { + // Do projection + { + List<pointConstraint> constraints(lambdas.size()); + pointField start(points); + searchableSurfacesQueries::findNearest + ( + geometry_, + surfaces_, + start, + scalarField(start.size(), magSqr(d)), + points, + constraints + ); + + // Reset start and end point + if (lambdas[0] < SMALL) + { + points[0] = startPt; + } + if (lambdas.last() > 1.0-SMALL) + { + points.last() = endPt; + } + + if (debugStr.valid()) + { + forAll(points, i) + { + debugStr().write(linePointRef(start[i], points[i])); + } + } + } + + // Calculate lambdas (normalised coordinate along edge) + scalarField projLambdas(points.size()); + { + projLambdas[0] = 0.0; + for (label i = 1; i < points.size(); i++) + { + projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]); + } + projLambdas /= projLambdas.last(); + } + linearInterpolationWeights interpolator(projLambdas); + + // Compare actual distances and move points (along straight line; + // not along surface) + vectorField residual(points.size(), vector::zero); + labelList indices; + scalarField weights; + for (label i = 1; i < points.size() - 1; i++) + { + interpolator.valueWeights(lambdas[i], indices, weights); + + point predicted = vector::zero; + forAll(indices, indexi) + { + predicted += weights[indexi]*points[indices[indexi]]; + } + residual[i] = predicted-points[i]; + } + + scalar scalarResidual = sum(mag(residual)); + + if (debug) + { + Pout<< "Iter:" << iter << " initialResidual:" << initialResidual + << " residual:" << scalarResidual << endl; + } + + if (scalarResidual < absTol*0.5*lambdas.size()) + { + break; + } + else if (iter == 0) + { + initialResidual = scalarResidual; + } + else if (scalarResidual/initialResidual < relTol) + { + break; + } + + + if (debugStr.valid()) + { + forAll(points, i) + { + const linePointRef ln(points[i], points[i]+residual[i]); + debugStr().write(ln); + } + } + + points += residual; + } + + return tpoints; +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H new file mode 100644 index 0000000000000000000000000000000000000000..27d7f7ffcac5bb860500e90cf912b86b5692c60a --- /dev/null +++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::projectCurveEdge + +Description + Defines the edge from the projection onto a surface (single surface) + or intersection of two surfaces. + +SourceFiles + projectCurveEdge.C + +\*---------------------------------------------------------------------------*/ + +#ifndef projectCurveEdge_H +#define projectCurveEdge_H + +#include "blockEdge.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class pointConstraint; + +/*---------------------------------------------------------------------------*\ + Class projectCurveEdge Declaration +\*---------------------------------------------------------------------------*/ + +class projectCurveEdge +: + public blockEdge +{ + // Private data + + const searchableSurfaces& geometry_; + + //- The indices of surfaces onto which the points are projected + labelList surfaces_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + projectCurveEdge(const projectCurveEdge&); + + //- Disallow default bitwise assignment + void operator=(const projectCurveEdge&); + + +public: + + //- Runtime type information + TypeName("projectCurve"); + + + // Constructors + + //- Construct from Istream setting pointsList + projectCurveEdge + ( + const dictionary& dict, + const label index, + const searchableSurfaces& geometry, + const pointField& points, + Istream& + ); + + + //- Destructor + virtual ~projectCurveEdge() + {} + + + // Member Functions + + //- Return the point positions corresponding to the curve parameters + // 0 <= lambda <= 1 + virtual point position(const scalar) const + { + NotImplemented; + return point::max; + } + + //- Return the point positions corresponding to the curve parameters + // 0 <= lambda <= 1 + virtual tmp<pointField> position(const scalarList&) const; + + //- Return the length of the curve + virtual scalar length() const + { + NotImplemented; + return 1; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C index 9f10f341918cad067e0705ab25a0abce66c4aa13..eefdae7af893a1759b9c0de1629bb75308af65be 100644 --- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C +++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C @@ -29,6 +29,7 @@ License #include "addToRunTimeSelectionTable.H" #include "pointConstraint.H" #include "OBJstream.H" +#include "linearInterpolationWeights.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -79,12 +80,14 @@ void Foam::projectEdge::findNearest Foam::projectEdge::projectEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ) : - blockEdge(points, is), + blockEdge(dict, index, points, is), geometry_(geometry) { wordList names(is); @@ -125,7 +128,20 @@ Foam::point Foam::projectEdge::position(const scalar lambda) const Foam::tmp<Foam::pointField> Foam::projectEdge::position(const scalarList& lambdas) const { - static label iter = 0; + // For debugging to tag the output + static label eIter = 0; + + autoPtr<OBJstream> debugStr; + if (debug) + { + debugStr.reset + ( + new OBJstream("projectEdge_" + Foam::name(eIter++) + ".obj") + ); + Info<< "Writing lines from straight-line start points" + << " to projected points to " << debugStr().name() << endl; + } + tmp<pointField> tpoints(new pointField(lambdas.size())); pointField& points = tpoints.ref(); @@ -141,17 +157,26 @@ Foam::projectEdge::position(const scalarList& lambdas) const } - for (label i = 0; i < 3; i++) + // Upper limit for number of iterations + const label maxIter = 10; + // Residual tolerance + const scalar relTol = 0.1; + const scalar absTol = 1e-4; + + scalar initialResidual = 0.0; + + for (label iter = 0; iter < maxIter; iter++) { // Do projection { List<pointConstraint> constraints(lambdas.size()); + pointField start(points); searchableSurfacesQueries::findNearest ( geometry_, surfaces_, - pointField(points), - scalarField(points.size(), magSqr(d)), + start, + scalarField(start.size(), magSqr(d)), points, constraints ); @@ -165,63 +190,77 @@ Foam::projectEdge::position(const scalarList& lambdas) const { points.last() = endPt; } + + if (debugStr.valid()) + { + forAll(points, i) + { + debugStr().write(linePointRef(start[i], points[i])); + } + } } - // Calculate distances - scalarField nearLength(points.size()); + // Calculate lambdas (normalised coordinate along edge) + scalarField projLambdas(points.size()); { - nearLength[0] = 0.0; - for(label i = 1; i < points.size(); i++) + projLambdas[0] = 0.0; + for (label i = 1; i < points.size(); i++) { - nearLength[i] = nearLength[i-1] + mag(points[i]-points[i-1]); + projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]); } + projLambdas /= projLambdas.last(); } + linearInterpolationWeights interpolator(projLambdas); // Compare actual distances and move points (along straight line; // not along surface) - for(label i = 1; i < points.size() - 1; i++) + vectorField residual(points.size(), vector::zero); + labelList indices; + scalarField weights; + for (label i = 1; i < points.size() - 1; i++) { - scalar nearDelta = mag(points[i]-points[i-1])/nearLength.last(); - scalar wantedDelta = lambdas[i]-lambdas[i-1]; + interpolator.valueWeights(lambdas[i], indices, weights); - vector v(points[i]-points[i-1]); - points[i] = points[i-1]+wantedDelta/nearDelta*v; + point predicted = vector::zero; + forAll(indices, indexi) + { + predicted += weights[indexi]*points[indices[indexi]]; + } + residual[i] = predicted-points[i]; } - } + scalar scalarResidual = sum(mag(residual)); - if (debug) - { - OBJstream str("projectEdge_" + Foam::name(iter++) + ".obj"); - Info<< "Writing lines from straight-line start points" - << " to projected points to " << str.name() << endl; - - pointField startPts(lambdas.size()); - forAll(lambdas, i) + if (debug) { - startPts[i] = startPt+lambdas[i]*d; + Pout<< "Iter:" << iter << " initialResidual:" << initialResidual + << " residual:" << scalarResidual << endl; } - pointField nearPts(lambdas.size()); - List<pointConstraint> nearConstraints(lambdas.size()); + if (scalarResidual < absTol*0.5*lambdas.size()) { - const scalar distSqr = magSqr(d); - searchableSurfacesQueries::findNearest - ( - geometry_, - surfaces_, - startPts, - scalarField(startPts.size(), distSqr), - nearPts, - nearConstraints - ); + break; + } + else if (iter == 0) + { + initialResidual = scalarResidual; + } + else if (scalarResidual/initialResidual < relTol) + { + break; } - forAll(startPts, i) + + if (debugStr.valid()) { - str.write(linePointRef(startPts[i], nearPts[i])); - str.write(linePointRef(nearPts[i], points[i])); + forAll(points, i) + { + const linePointRef ln(points[i], points[i]+residual[i]); + debugStr().write(ln); + } } + + points += residual; } return tpoints; diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H index cf8baebe75edf3a8d824d3e099e2305c2f558165..320f658eb57411b89e1ed39deb8ecf1163960534 100644 --- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H +++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H @@ -84,6 +84,8 @@ public: //- Construct from Istream setting pointsList projectEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C index 5a6be74f4e8dd94b6b87d1fd777825f72eb069de..0e70e255c82b2b74a44e2093ffe46712387fe3b9 100644 --- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C +++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C @@ -62,12 +62,14 @@ Foam::blockEdges::splineEdge::splineEdge Foam::blockEdges::splineEdge::splineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField& points, Istream& is ) : - blockEdge(points, is), + blockEdge(dict, index, points, is), CatmullRomSpline(appendEndPoints(points, start_, end_, pointField(is))) { token t(is); diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H index 4080b606c82a7b72698632567851e96c5ecc5da5..e1457a2d87070597293918e400b6f26bca495917 100644 --- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H +++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H @@ -83,6 +83,8 @@ public: //- Construct from Istream, setting pointsList splineEdge ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, const pointField&, Istream& diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C index 388c9fcc5dc9402c8e11d4fc9a8135584d7875d6..4ee92176c76fbac769a8e7c29100d6aa8aa58527 100644 --- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C +++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "blockFace.H" +#include "blockDescriptor.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -42,9 +43,21 @@ Foam::blockFace::blockFace(const face& vertices) {} -Foam::blockFace::blockFace(Istream& is) +Foam::blockFace::blockFace +( + const dictionary& dict, + const label index, + Istream& is +) : - vertices_(is) + vertices_ + ( + blockDescriptor::read<label> + ( + is, + dict.subOrEmptyDict("namedVertices") + ) + ) {} @@ -57,6 +70,8 @@ Foam::autoPtr<Foam::blockFace> Foam::blockFace::clone() const Foam::autoPtr<Foam::blockFace> Foam::blockFace::New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ) @@ -81,7 +96,36 @@ Foam::autoPtr<Foam::blockFace> Foam::blockFace::New << abort(FatalError); } - return autoPtr<blockFace>(cstrIter()(geometry, is)); + return autoPtr<blockFace>(cstrIter()(dict, index, geometry, is)); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::blockFace::write(Ostream& os, const dictionary& d) const +{ + const dictionary* varDictPtr = d.subDictPtr("namedVertices"); + if (varDictPtr) + { + const dictionary& varDict = *varDictPtr; + + // Write size and start delimiter + os << vertices_.size() << token::BEGIN_LIST; + + // Write contents + forAll(vertices_, i) + { + if (i > 0) os << token::SPACE; + blockDescriptor::write(os, vertices_[i], varDict); + } + + // Write end delimiter + os << token::END_LIST; + } + else + { + os << vertices_ << endl; + } } diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H index fee6771d79bb6b6a2dd08cef7626aa180a675bb6..f890aa3ae79bf79a08335a35f648a770184c6d26 100644 --- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H +++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H @@ -77,10 +77,12 @@ public: blockFace, Istream, ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ), - (geometry, is) + (dict, index, geometry, is) ); @@ -90,7 +92,12 @@ public: blockFace(const face& vertices); //- Construct from Istream - blockFace(Istream&); + blockFace + ( + const dictionary& dict, + const label index, + Istream& + ); //- Clone function virtual autoPtr<blockFace> clone() const; @@ -98,6 +105,8 @@ public: //- New function which constructs and returns pointer to a blockFace static autoPtr<blockFace> New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& ); @@ -106,18 +115,22 @@ public: // PtrLists of blockFace class iNew { + const dictionary& dict_; const searchableSurfaces& geometry_; + mutable label index_; public: - iNew(const searchableSurfaces& geometry) + iNew(const dictionary& dict, const searchableSurfaces& geometry) : - geometry_(geometry) + dict_(dict), + geometry_(geometry), + index_(0) {} autoPtr<blockFace> operator()(Istream& is) const { - return blockFace::New(geometry_, is); + return blockFace::New(dict_, index_++, geometry_, is); } }; @@ -145,6 +158,9 @@ public: pointField& points ) const = 0; + //- Write face with variable backsubstitution + void write(Ostream&, const dictionary&) const; + // Ostream operator diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C index 2cf9c9742769e71c72f30e3c6cf9fd86b2e60cc0..13ca420a1c7480b599a682ed21b410b463300bff 100644 --- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C +++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C @@ -27,6 +27,8 @@ License #include "unitConversion.H" #include "addToRunTimeSelectionTable.H" #include "blockDescriptor.H" +#include "OBJstream.H" +#include "linearInterpolationWeights.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -66,15 +68,74 @@ const Foam::searchableSurface& Foam::blockFaces::projectFace::lookupSurface } +Foam::label Foam::blockFaces::projectFace::index +( + const labelPair& n, + const labelPair& coord +) const +{ + return coord.first()+coord.second()*n.first(); +} + + +void Foam::blockFaces::projectFace::calcLambdas +( + const labelPair& n, + const pointField& points, + scalarField& lambdaI, + scalarField& lambdaJ +) const +{ + lambdaI.setSize(points.size()); + lambdaI = 0.0; + lambdaJ.setSize(points.size()); + lambdaJ = 0.0; + + for (label i = 1; i < n.first(); i++) + { + for (label j = 1; j < n.second(); j++) + { + label ij = index(n, labelPair(i, j)); + label iMin1j = index(n, labelPair(i-1, j)); + lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij]-points[iMin1j]); + + label ijMin1 = index(n, labelPair(i, j-1)); + lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij]-points[ijMin1]); + } + } + + for (label i = 1; i < n.first(); i++) + { + label ijLast = index(n, labelPair(i, n.second()-1)); + for (label j = 1; j < n.second(); j++) + { + label ij = index(n, labelPair(i, j)); + lambdaJ[ij] /= lambdaJ[ijLast]; + } + } + for (label j = 1; j < n.second(); j++) + { + label iLastj = index(n, labelPair(n.first()-1, j)); + for (label i = 1; i < n.first(); i++) + { + label ij = index(n, labelPair(i, j)); + lambdaI[ij] /= lambdaI[iLastj]; + } + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::blockFaces::projectFace::projectFace ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ) : - blockFace(is), + blockFace(dict, index, is), surface_(lookupSurface(geometry, is)) {} @@ -88,20 +149,215 @@ void Foam::blockFaces::projectFace::project pointField& points ) const { - List<pointIndexHit> hits; - scalarField nearestDistSqr - ( - points.size(), - magSqr(points[0] - points[points.size()-1]) - ); - surface_.findNearest(points, nearestDistSqr, hits); - - forAll(hits, i) + // For debugging to tag the output + static label fIter = 0; + + autoPtr<OBJstream> debugStr; + if (debug) { - if (hits[i].hit()) + debugStr.reset + ( + new OBJstream("projectFace_" + Foam::name(fIter++) + ".obj") + ); + Info<< "Face:" << blockFacei << " on block:" << desc.blockShape() + << " with verts:" << desc.vertices() + << " writing lines from start points" + << " to projected points to " << debugStr().name() << endl; + } + + + // Find out range of vertices in face + labelPair n(-1, -1); + switch (blockFacei) + { + case 0: + case 1: + { + n.first() = desc.density()[1]+1; + n.second() = desc.density()[2]+1; + } + break; + + case 2: + case 3: + { + n.first() = desc.density()[0]+1; + n.second() = desc.density()[2]+1; + } + break; + + case 4: + case 5: { - points[i] = hits[i].hitPoint(); + n.first() = desc.density()[0]+1; + n.second() = desc.density()[1]+1; } + break; + } + + + // Calculate initial normalised edge lengths (= u,v coordinates) + scalarField lambdaI(points.size(), 0.0); + scalarField lambdaJ(points.size(), 0.0); + calcLambdas(n, points, lambdaI, lambdaJ); + + + // Upper limit for number of iterations + const label maxIter = 10; + // Residual tolerance + const scalar relTol = 0.1; + + scalar initialResidual = 0.0; + scalar iResidual = 0.0; + scalar jResidual = 0.0; + + for (label iter = 0; iter < maxIter; iter++) + { + // Do projection + { + List<pointIndexHit> hits; + scalarField nearestDistSqr + ( + points.size(), + magSqr(points[0] - points[points.size()-1]) + ); + surface_.findNearest(points, nearestDistSqr, hits); + + forAll(hits, i) + { + if (hits[i].hit()) + { + const point& hitPt = hits[i].hitPoint(); + if (debugStr.valid()) + { + debugStr().write(linePointRef(points[i], hitPt)); + } + points[i] = hitPt; + } + } + } + + if (debug) + { + Pout<< "Iter:" << iter << " initialResidual:" << initialResidual + << " iResidual+jResidual:" << iResidual+jResidual << endl; + } + + + if (iter > 0 && (iResidual+jResidual)/initialResidual < relTol) + { + break; + } + + + // Predict along i + vectorField residual(points.size(), vector::zero); + + // Work arrays for interpolation + labelList indices; + scalarField weights; + for (label j = 1; j < n.second()-1; j++) + { + // Calculate actual lamdba along constant j line + scalarField projLambdas(n.first()); + projLambdas[0] = 0.0; + for (label i = 1; i < n.first(); i++) + { + label ij = index(n, labelPair(i, j)); + label iMin1j = index(n, labelPair(i-1, j)); + projLambdas[i] = + projLambdas[i-1] + +mag(points[ij]-points[iMin1j]); + } + projLambdas /= projLambdas.last(); + + linearInterpolationWeights interpolator(projLambdas); + + for (label i = 1; i < n.first()-1; i++) + { + label ij = index(n, labelPair(i, j)); + + interpolator.valueWeights(lambdaI[ij], indices, weights); + + point predicted = vector::zero; + forAll(indices, indexi) + { + label ptIndex = index(n, labelPair(indices[indexi], j)); + predicted += weights[indexi]*points[ptIndex]; + } + residual[ij] = predicted-points[ij]; + } + } + + if (debugStr.valid()) + { + forAll(points, i) + { + const linePointRef ln(points[i], points[i]+residual[i]); + debugStr().write(ln); + } + } + + iResidual = sum(mag(residual)); + + // Update points before doing j. Note: is this needed? Complicates + // residual checking. + points += residual; + + + // Predict along j + residual = vector::zero; + for (label i = 1; i < n.first()-1; i++) + { + // Calculate actual lamdba along constant i line + scalarField projLambdas(n.second()); + projLambdas[0] = 0.0; + for (label j = 1; j < n.second(); j++) + { + label ij = index(n, labelPair(i, j)); + label ijMin1 = index(n, labelPair(i, j-1)); + projLambdas[j] = + projLambdas[j-1] + +mag(points[ij]-points[ijMin1]); + } + + projLambdas /= projLambdas.last(); + + linearInterpolationWeights interpolator(projLambdas); + + for (label j = 1; j < n.second()-1; j++) + { + label ij = index(n, labelPair(i, j)); + + interpolator.valueWeights(lambdaJ[ij], indices, weights); + + point predicted = vector::zero; + forAll(indices, indexi) + { + label ptIndex = index(n, labelPair(i, indices[indexi])); + predicted += weights[indexi]*points[ptIndex]; + } + residual[ij] = predicted-points[ij]; + } + } + + if (debugStr.valid()) + { + forAll(points, i) + { + const linePointRef ln(points[i], points[i]+residual[i]); + debugStr().write(ln); + } + } + + jResidual = sum(mag(residual)); + + if (iter == 0) + { + initialResidual = iResidual + jResidual; + } + + points += residual; } } diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H index c1b56bbcdba83e101bc1ae23be684b33d011e5dd..635318b05c27501c95bdb281a87a94e1f9f33409 100644 --- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H +++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H @@ -67,6 +67,22 @@ class projectFace Istream& is ) const; + //- Convert i,j to single index + label index + ( + const labelPair& n, + const labelPair& coord + ) const; + + //- Calulate lambdas (but unnormalised) + void calcLambdas + ( + const labelPair& n, + const pointField& points, + scalarField& lambdaI, + scalarField& lambdaJ + ) const; + //- Disallow default bitwise copy construct projectFace(const projectFace&); @@ -85,6 +101,8 @@ public: //- Construct from Istream setting pointsList projectFace ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& ); diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C index 366288e5ae72ef566d50681c0262a60650a2844f..4a08adbfbee76cf8209a943e27aa176f272ebcc3 100644 --- a/src/mesh/blockMesh/blockMesh/blockMesh.C +++ b/src/mesh/blockMesh/blockMesh/blockMesh.C @@ -59,7 +59,7 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName) blockVertices_ ( dict.lookup("vertices"), - blockVertex::iNew(geometry_) + blockVertex::iNew(dict, geometry_) ), vertices_(Foam::vertices(blockVertices_)), topologyPtr_(createTopology(dict, regionName)) diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H index fcd3698c34027ee3f08cf42be98ce2751a5d25af..093972d12d02afd29113f5aad51d60241b40c04e 100644 --- a/src/mesh/blockMesh/blockMesh/blockMesh.H +++ b/src/mesh/blockMesh/blockMesh/blockMesh.H @@ -139,7 +139,7 @@ class blockMesh polyMesh* createTopology(const IOdictionary&, const word& regionName); - void check(const polyMesh&) const; + void check(const polyMesh&, const dictionary&) const; //- Determine the merge info and the final number of cells/points void calcMergeInfo(); diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C index 10bd9e3116df63d26c852f0972749ffa94da3752..fc0c08b8623443bfcd277132f852d8755a958a00 100644 --- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C +++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C @@ -27,7 +27,7 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::blockMesh::check(const polyMesh& bm) const +void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const { Info<< nl << "Check topology" << endl; @@ -40,8 +40,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const { if (edges_[cei].compare(edges_[cej]) != 0) { - Info<< " Curved edge " << edges_[cej] - << " is a duplicate of curved edge " << edges_[cei] + Info<< " Curved edge "; + edges_[cej].write(Info, dict); + Info<< " is a duplicate of curved edge " << edges_[cei] << endl; ok = false; break; @@ -74,8 +75,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const if (!found) { - Info<< " Curved edge " << edges_[cei] - << " does not correspond to a block edge." + Info<< " Curved edge "; + edges_[cei].write(Info, dict); + Info<< " does not correspond to a block edge." << endl; ok = false; } @@ -90,9 +92,11 @@ void Foam::blockMesh::check(const polyMesh& bm) const { if (faces_[cfi].compare(faces_[cfj]) != 0) { - Info<< " Curved face " << faces_[cfj] - << " is a duplicate of curved face " << faces_[cfi] - << endl; + Info<< " Curved face "; + faces_[cfj].write(Info, dict); + Info<< " is a duplicate of curved face "; + faces_[cfi].write(Info, dict); + Info<< endl; ok = false; break; } @@ -112,9 +116,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const if (!found) { - Info<< " Curved face " << faces_[cfi] - << " does not correspond to a block face." - << endl; + Info<< " Curved face "; + faces_[cfi].write(Info, dict); + Info<< " does not correspond to a block face." << endl; ok = false; } } diff --git a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C index ab998442051ad93e6b1bf0bc884573f33f857f43..a4c9acead3f9477f02e30cfe532f6c89c8518c9e 100644 --- a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C +++ b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C @@ -110,6 +110,11 @@ void Foam::blockMesh::readPatches wordList& nbrPatchNames ) { + // Collect all variables + dictionary varDict(meshDescription.subOrEmptyDict("namedVertices")); + varDict.merge(meshDescription.subOrEmptyDict("namedBlocks")); + + ITstream& patchStream(meshDescription.lookup("patches")); // Read number of patches in mesh @@ -160,7 +165,11 @@ void Foam::blockMesh::readPatches >> patchNames[nPatches]; // Read patch faces - patchStream >> tmpBlocksPatches[nPatches]; + tmpBlocksPatches[nPatches] = blockDescriptor::read<face> + ( + patchStream, + varDict + ); // Check for multiple patches @@ -258,6 +267,11 @@ void Foam::blockMesh::readBoundary PtrList<dictionary>& patchDicts ) { + // Collect all variables + dictionary varDict(meshDescription.subOrEmptyDict("namedVertices")); + varDict.merge(meshDescription.subOrEmptyDict("namedBlocks")); + + // Read like boundary file const PtrList<entry> patchesInfo ( @@ -286,7 +300,12 @@ void Foam::blockMesh::readBoundary patchDicts.set(patchi, new dictionary(patchInfo.dict())); // Read block faces - patchDicts[patchi].lookup("faces") >> tmpBlocksPatches[patchi]; + tmpBlocksPatches[patchi] = blockDescriptor::read<face> + ( + patchDicts[patchi].lookup("faces"), + varDict + ); + checkPatchLabels ( @@ -353,7 +372,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology blockEdgeList edges ( meshDescription.lookup("edges"), - blockEdge::iNew(geometry_, vertices_) + blockEdge::iNew(meshDescription, geometry_, vertices_) ); edges_.transfer(edges); @@ -375,7 +394,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology blockFaceList faces ( meshDescription.lookup("faces"), - blockFace::iNew(geometry_) + blockFace::iNew(meshDescription, geometry_) ); faces_.transfer(faces); @@ -395,7 +414,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology blockList blocks ( meshDescription.lookup("blocks"), - block::iNew(vertices_, edges_, faces_) + block::iNew(meshDescription, vertices_, edges_, faces_) ); transfer(blocks); @@ -545,7 +564,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology ); } - check(*blockMeshPtr); + check(*blockMeshPtr, meshDescription); return blockMeshPtr; } diff --git a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C index fe8ab41dc590827f32ef1f0a22e6b9586c83871c..6bdf9c2adee47f3297b75bb891de0c90f0ac2a5a 100644 --- a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C +++ b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C @@ -50,6 +50,8 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::clone() const Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ) @@ -61,14 +63,14 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New token firstToken(is); - if (firstToken.pToken() == token::BEGIN_LIST) + if (firstToken.isPunctuation() && firstToken.pToken() == token::BEGIN_LIST) { // Putback the opening bracket is.putBack(firstToken); return autoPtr<blockVertex> ( - new blockVertices::pointVertex(geometry, is) + new blockVertices::pointVertex(dict, index, geometry, is) ); } else if (firstToken.isWord()) @@ -88,7 +90,7 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New << abort(FatalError); } - return autoPtr<blockVertex>(cstrIter()(geometry, is)); + return autoPtr<blockVertex>(cstrIter()(dict, index, geometry, is)); } else { diff --git a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H index 7d5bb9fbd635c16ca9a6da3774578f8b7640fdab..5c6c8006a1b27366a4d283cac1cc41464e57a02e 100644 --- a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H +++ b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H @@ -63,10 +63,12 @@ public: blockVertex, Istream, ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ), - (geometry, is) + (dict, index, geometry, is) ); @@ -81,6 +83,8 @@ public: //- New function which constructs and returns pointer to a blockVertex static autoPtr<blockVertex> New ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& ); @@ -89,18 +93,22 @@ public: // PtrLists of blockVertex class iNew { + const dictionary& dict_; const searchableSurfaces& geometry_; + mutable label index_; public: - iNew(const searchableSurfaces& geometry) + iNew(const dictionary& dict, const searchableSurfaces& geometry) : - geometry_(geometry) + dict_(dict), + geometry_(geometry), + index_(0) {} autoPtr<blockVertex> operator()(Istream& is) const { - return blockVertex::New(geometry_, is); + return blockVertex::New(dict_, index_++, geometry_, is); } }; diff --git a/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C new file mode 100644 index 0000000000000000000000000000000000000000..b4392fac2689cf2c824f69ff1ef62bd016c08cef --- /dev/null +++ b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C @@ -0,0 +1,81 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "namedVertex.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace blockVertices +{ + defineTypeNameAndDebug(namedVertex, 0); + addToRunTimeSelectionTable(blockVertex, namedVertex, Istream); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::blockVertices::namedVertex::namedVertex +( + const dictionary& dict, + const label index, + const searchableSurfaces& geometry, + Istream& is +) +: + name_(is), + vertexPtr_(blockVertex::New(dict, index, geometry, is)) +{ + //Info<< "Vertex " << name_ << " at " << vertexPtr_().operator point() + // << " has index " << index << endl; + + dictionary& d = const_cast<dictionary&>(dict); + + const dictionary* varDictPtr = d.subDictPtr("namedVertices"); + if (varDictPtr) + { + const_cast<dictionary&>(*varDictPtr).add(name_, index); + } + else + { + dictionary varDict; + varDict.add(name_, index); + d.add("namedVertices", varDict); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::blockVertices::namedVertex::operator point() const +{ + return vertexPtr_().operator point(); +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H new file mode 100644 index 0000000000000000000000000000000000000000..117107ad5f5fe12e49a2f449b1113e0480cffce9 --- /dev/null +++ b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::blockVertices::namedVertex + +Description + Gives name to a vertex. + +SourceFiles + namedVertex.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockVertices_namedVertex_H +#define blockVertices_namedVertex_H + +#include "blockVertex.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace blockVertices +{ + +/*---------------------------------------------------------------------------*\ + Class namedVertex Declaration +\*---------------------------------------------------------------------------*/ + +class namedVertex +: + public blockVertex +{ +protected: + + // Protected member data + + //- The dictionary variable name for the vertex number + word name_; + + //- The vertex location + autoPtr<blockVertex> vertexPtr_; + + +public: + + //- Runtime type information + TypeName("named"); + + + // Constructors + + //- Construct from Istream setting pointsList + namedVertex + ( + const dictionary&, + const label index, + const searchableSurfaces& geometry, + Istream& + ); + + + //- Destructor + virtual ~namedVertex() + {} + + + // Member Functions + + virtual operator point() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace blockVertices +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C index 4b3120d98b2a3bea482f4366650d0ec5293ad573..4af95df1637f5aec1a8fe61a4a96a055dda5deb2 100644 --- a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C +++ b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C @@ -42,6 +42,8 @@ namespace blockVertices Foam::blockVertices::pointVertex::pointVertex ( + const dictionary&, + const label index, const searchableSurfaces& geometry, Istream& is ) diff --git a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H index 1179c66bd5d8e2da6c70d9f7cb68276f1f0e6d21..929f92a25c201acc8caebb0ebd0da562923100a5 100644 --- a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H +++ b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H @@ -70,6 +70,8 @@ public: //- Construct from Istream setting pointsList pointVertex ( + const dictionary&, + const label index, const searchableSurfaces& geometry, Istream& ); diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C index 04cff7d250343172f36a0b60ea58c2135aed0cba..f5d23fc50c9f1c21cc92a9678c4970c742ce3a28 100644 --- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C +++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C @@ -45,11 +45,13 @@ namespace blockVertices Foam::blockVertices::projectVertex::projectVertex ( + const dictionary& dict, + const label index, const searchableSurfaces& geometry, Istream& is ) : - pointVertex(geometry, is), + pointVertex(dict, index, geometry, is), geometry_(geometry) { wordList names(is); @@ -79,8 +81,11 @@ Foam::blockVertices::projectVertex::operator point() const // Note: how far do we need to search? Probably not further than - // span of surfaces themselves. + // span of surfaces themselves. Make sure to limit in case + // of e.g. searchablePlane which has infinite bb. boundBox bb(searchableSurfacesQueries::bounds(geometry_, surfaces_)); + bb.min() = max(bb.min(), point(-GREAT, -GREAT, -GREAT)); + bb.max() = min(bb.max(), point(GREAT, GREAT, GREAT)); searchableSurfacesQueries::findNearest ( diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H index 929b6a68ebe646fb6136e60cf3b520d5b517af42..d47dce6e8cbd89403986fce2f31f9b2c266bda3f 100644 --- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H +++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H @@ -81,6 +81,8 @@ public: //- Construct from Istream setting pointsList projectVertex ( + const dictionary&, + const label index, const searchableSurfaces& geometry, Istream& ); diff --git a/src/mesh/blockMesh/block/block.C b/src/mesh/blockMesh/blocks/block/block.C similarity index 61% rename from src/mesh/blockMesh/block/block.C rename to src/mesh/blockMesh/blocks/block/block.C index 1722b0095ee0f51b8df67362ea10459fff768e9a..5d75aa5574f40903e41c973e6cd79b3623252eb2 100644 --- a/src/mesh/blockMesh/block/block.C +++ b/src/mesh/blockMesh/blocks/block/block.C @@ -25,17 +25,27 @@ License #include "block.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(block, 0); + defineRunTimeSelectionTable(block, Istream); +} + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::block::block ( + const dictionary& dict, + const label index, const pointField& vertices, const blockEdgeList& edges, const blockFaceList& faces, Istream& is ) : - blockDescriptor(vertices, edges, faces, is) + blockDescriptor(dict, index, vertices, edges, faces, is) { createPoints(); createBoundary(); @@ -51,6 +61,49 @@ Foam::block::block(const blockDescriptor& blockDesc) } +Foam::autoPtr<Foam::block> Foam::block::New +( + const dictionary& dict, + const label index, + const pointField& points, + const blockEdgeList& edges, + const blockFaceList& faces, + Istream& is +) +{ + if (debug) + { + InfoInFunction << "Constructing block" << endl; + } + + const word blockOrCellShapeType(is); + + IstreamConstructorTable::iterator cstrIter = + IstreamConstructorTablePtr_->find(blockOrCellShapeType); + + if (cstrIter == IstreamConstructorTablePtr_->end()) + { + is.putBack(blockOrCellShapeType); + return autoPtr<block>(new block(dict, index, points, edges, faces, is)); + } + else + { + return autoPtr<block> + ( + cstrIter() + ( + dict, + index, + points, + edges, + faces, + is + ) + ); + } +} + + // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const block& b) diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/blocks/block/block.H similarity index 78% rename from src/mesh/blockMesh/block/block.H rename to src/mesh/blockMesh/blocks/block/block.H index 7b55c5f59f1afb29765d7969b8b6e5170b0e9445..fb3012016d55d8b9ee1f5ad9789158c0ea4d4069 100644 --- a/src/mesh/blockMesh/block/block.H +++ b/src/mesh/blockMesh/blocks/block/block.H @@ -91,11 +91,36 @@ class block public: + //- Runtime type information + TypeName("block"); + + + // Declare run-time constructor selection tables + + declareRunTimeSelectionTable + ( + autoPtr, + block, + Istream, + ( + const dictionary& dict, + const label index, + const pointField& vertices, + const blockEdgeList& edges, + const blockFaceList& faces, + Istream& is + ), + (dict, index, vertices, edges, faces, is) + ); + + // Constructors //- Construct from components with Istream block ( + const dictionary& dict, + const label index, const pointField& vertices, const blockEdgeList& edges, const blockFaceList& faces, @@ -112,35 +137,56 @@ public: return autoPtr<block>(nullptr); } + //- New function which constructs and returns pointer to a block + static autoPtr<block> New + ( + const dictionary& dict, + const label index, + const pointField& points, + const blockEdgeList& edges, + const blockFaceList& faces, + Istream& + ); + //- Class used for the read-construction of // PtrLists of blocks class iNew { + const dictionary& dict_; const pointField& points_; const blockEdgeList& edges_; const blockFaceList& faces_; + mutable label index_; public: iNew ( + const dictionary& dict, const pointField& points, const blockEdgeList& edges, const blockFaceList& faces ) : + dict_(dict), points_(points), edges_(edges), - faces_(faces) + faces_(faces), + index_(0) {} autoPtr<block> operator()(Istream& is) const { - return autoPtr<block>(new block(points_, edges_, faces_, is)); + return block::New(dict_, index_++, points_, edges_, faces_, is); } }; + //- Destructor + virtual ~block() + {} + + // Member Functions // Access diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/blocks/block/blockCreate.C similarity index 100% rename from src/mesh/blockMesh/block/blockCreate.C rename to src/mesh/blockMesh/blocks/block/blockCreate.C diff --git a/src/mesh/blockMesh/block/blockI.H b/src/mesh/blockMesh/blocks/block/blockI.H similarity index 100% rename from src/mesh/blockMesh/block/blockI.H rename to src/mesh/blockMesh/blocks/block/blockI.H diff --git a/src/mesh/blockMesh/block/blockList.H b/src/mesh/blockMesh/blocks/block/blockList.H similarity index 95% rename from src/mesh/blockMesh/block/blockList.H rename to src/mesh/blockMesh/blocks/block/blockList.H index bfca6c9c60848ccaad05a9644190b634946b7ac4..80b8a328af597b678718eb4c97c456d05a2a860e 100644 --- a/src/mesh/blockMesh/block/blockList.H +++ b/src/mesh/blockMesh/blocks/block/blockList.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -53,4 +53,3 @@ typedef PtrList<block> blockList; #endif // ************************************************************************* // - diff --git a/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C new file mode 100644 index 0000000000000000000000000000000000000000..9408e4fd3e82eb20f5cb1a19b123dfcb1d1ee763 --- /dev/null +++ b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "namedBlock.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace blocks +{ + defineTypeNameAndDebug(namedBlock, 0); + addToRunTimeSelectionTable(block, namedBlock, Istream); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::blocks::namedBlock::namedBlock +( + const dictionary& dict, + const label index, + const pointField& vertices, + const blockEdgeList& edges, + const blockFaceList& faces, + Istream& is +) +: + word(is), + block(dict, index, vertices, edges, faces, is) +{ + //Info<< "Block " << static_cast<word&>(*this) + // << " has index " << index << endl; + + dictionary& d = const_cast<dictionary&>(dict); + + const dictionary* varDictPtr = d.subDictPtr("namedBlocks"); + if (varDictPtr) + { + const_cast<dictionary&>(*varDictPtr).add(*this, index); + } + else + { + dictionary varDict; + varDict.add(*this, index); + d.add("namedBlocks", varDict); + } +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H new file mode 100644 index 0000000000000000000000000000000000000000..98968f60b0659e6b125dd75b730d1d7d9f242c54 --- /dev/null +++ b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::blocks::namedBlock + +Description + Gives name to a block. + +SourceFiles + namedBlock.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blocks_namedBlock_H +#define blocks_namedBlock_H + +#include "block.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace blocks +{ + +/*---------------------------------------------------------------------------*\ + Class namedBlock Declaration +\*---------------------------------------------------------------------------*/ + +class namedBlock +: + public word, + public block +{ +public: + + //- Runtime type information + TypeName("named"); + + + // Constructors + + //- Construct from Istream setting pointsList + namedBlock + ( + const dictionary& dict, + const label index, + const pointField& vertices, + const blockEdgeList& edges, + const blockFaceList& faces, + Istream& is + ); + + + //- Destructor + virtual ~namedBlock() + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace blocks +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/searchableCurve/searchableCurve.C b/src/mesh/blockMesh/searchableCurve/searchableCurve.C new file mode 100644 index 0000000000000000000000000000000000000000..b8cc234d0e8c72de7107b83c5326bfc4bc466329 --- /dev/null +++ b/src/mesh/blockMesh/searchableCurve/searchableCurve.C @@ -0,0 +1,432 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "searchableCurve.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "edgeMesh.H" +#include "indexedOctree.H" +#include "treeDataEdge.H" +#include "linearInterpolationWeights.H" +#include "quaternion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(searchableCurve, 0); + addToRunTimeSelectionTable(searchableSurface, searchableCurve, dict); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::searchableCurve::searchableCurve +( + const IOobject& io, + const dictionary& dict +) +: + searchableSurface(io), + eMeshPtr_ + ( + edgeMesh::New + ( + IOobject + ( + dict.lookup("file"), // name + io.time().constant(), // instance + "triSurface", // local + io.time(), // registry + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ).objectPath() + ) + ), + radius_(readScalar(dict.lookup("radius"))) +{ + const edgeMesh& eMesh = eMeshPtr_(); + + const pointField& points = eMesh.points(); + const edgeList& edges = eMesh.edges(); + bounds() = boundBox(points, false); + + vector halfSpan(0.5*bounds().span()); + point ctr(bounds().midpoint()); + + bounds().min() = ctr - mag(halfSpan)*vector(1, 1, 1); + bounds().max() = ctr + mag(halfSpan)*vector(1, 1, 1); + + // Calculate bb of all points + treeBoundBox bb(bounds()); + + // Slightly extended bb. Slightly off-centred just so on symmetric + // geometry there are less face/edge aligned items. + bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + + edgeTree_.reset + ( + new indexedOctree<treeDataEdge> + ( + treeDataEdge + ( + false, // do not cache bb + edges, + points, + identity(edges.size()) + ), + bb, // overall search domain + 8, // maxLevel + 10, // leafsize + 3.0 // duplicity + ) + ); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::searchableCurve::~searchableCurve() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::wordList& Foam::searchableCurve::regions() const +{ + if (regions_.empty()) + { + regions_.setSize(1); + regions_[0] = "region0"; + } + return regions_; +} + + +Foam::label Foam::searchableCurve::size() const +{ + return eMeshPtr_().points().size(); +} + + +Foam::tmp<Foam::pointField> Foam::searchableCurve::coordinates() const +{ + return eMeshPtr_().points(); +} + + +void Foam::searchableCurve::boundingSpheres +( + pointField& centres, + scalarField& radiusSqr +) const +{ + centres = eMeshPtr_().points(); + radiusSqr.setSize(centres.size()); + radiusSqr = Foam::sqr(radius_); + // Add a bit to make sure all points are tested inside + radiusSqr += Foam::sqr(SMALL); +} + + +void Foam::searchableCurve::findNearest +( + const pointField& samples, + const scalarField& nearestDistSqr, + List<pointIndexHit>& info +) const +{ + const indexedOctree<treeDataEdge>& tree = edgeTree_(); + + info.setSize(samples.size()); + + forAll(samples, i) + { + info[i] = tree.findNearest(samples[i], nearestDistSqr[i]); + + if (info[i].hit()) + { + vector d(samples[i]-info[i].hitPoint()); + info[i].setPoint(info[i].hitPoint() + d/mag(d)*radius_); + } + } +} + + +void Foam::searchableCurve::findNearest +( + const point& start, + const point& end, + const scalarField& rawLambdas, + const scalarField& nearestDistSqr, + List<pointIndexHit>& info +) const +{ + const edgeMesh& mesh = eMeshPtr_(); + const indexedOctree<treeDataEdge>& tree = edgeTree_(); + const edgeList& edges = mesh.edges(); + const pointField& points = mesh.points(); + const labelListList& pointEdges = mesh.pointEdges(); + + const scalar maxDistSqr(Foam::magSqr(bounds().span())); + + // Normalise lambdas + const scalarField lambdas + ( + (rawLambdas-rawLambdas[0]) + /(rawLambdas.last()-rawLambdas[0]) + ); + + + // Nearest point on curve and local axis direction + pointField curvePoints(lambdas.size()); + vectorField axialVecs(lambdas.size()); + + const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr); + curvePoints[0] = startInfo.hitPoint(); + axialVecs[0] = edges[startInfo.index()].vec(points); + + const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr); + curvePoints.last() = endInfo.hitPoint(); + axialVecs.last() = edges[endInfo.index()].vec(points); + + + + scalarField curveLambdas(points.size(), -1.0); + + { + scalar endDistance = -1.0; + + // Determine edge lengths walking from start to end. + + const point& start = curvePoints[0]; + const point& end = curvePoints.last(); + + label edgei = startInfo.index(); + const edge& startE = edges[edgei]; + + label pointi = startE[0]; + if ((startE.vec(points)&(end-start)) > 0) + { + pointi = startE[1]; + } + + curveLambdas[pointi] = mag(points[pointi]-curvePoints[0]); + label otherPointi = startE.otherVertex(pointi); + curveLambdas[otherPointi] = -mag(points[otherPointi]-curvePoints[0]); + + //Pout<< "for point:" << points[pointi] << " have distance " + // << curveLambdas[pointi] << endl; + + + while (true) + { + const labelList& pEdges = pointEdges[pointi]; + if (pEdges.size() == 1) + { + break; + } + else if (pEdges.size() != 2) + { + FatalErrorInFunction << "Curve " << name() + << " is not a single path. This is not supported" + << exit(FatalError); + break; + } + + label oldEdgei = edgei; + if (pEdges[0] == oldEdgei) + { + edgei = pEdges[1]; + } + else + { + edgei = pEdges[0]; + } + + if (edgei == endInfo.index()) + { + endDistance = curveLambdas[pointi] + mag(end-points[pointi]); + + //Pout<< "Found end edge:" << edges[edgei].centre(points) + // << " endPt:" << end + // << " point before:" << points[pointi] + // << " accumulated length:" << endDistance << endl; + } + + + label oldPointi = pointi; + pointi = edges[edgei].otherVertex(oldPointi); + + if (curveLambdas[pointi] >= 0) + { + break; + } + + curveLambdas[pointi] = + curveLambdas[oldPointi] + edges[edgei].mag(points); + } + + // Normalise curveLambdas + forAll(curveLambdas, i) + { + if (curveLambdas[i] >= 0) + { + curveLambdas[i] /= endDistance; + } + } + } + + + + // Interpolation engine + linearInterpolationWeights interpolator(curveLambdas); + + // Find wanted location along curve + labelList indices; + scalarField weights; + for (label i = 1; i < curvePoints.size()-1; i++) + { + interpolator.valueWeights(lambdas[i], indices, weights); + + if (indices.size() == 1) + { + // On outside of curve. Choose one of the connected edges. + label pointi = indices[0]; + const point& p0 = points[pointi]; + label edge0 = pointEdges[pointi][0]; + const edge& e0 = edges[edge0]; + axialVecs[i] = e0.vec(points); + curvePoints[i] = weights[0]*p0; + } + else if (indices.size() == 2) + { + const point& p0 = points[indices[0]]; + const point& p1 = points[indices[1]]; + axialVecs[i] = p1-p0; + curvePoints[i] = weights[0]*p0+weights[1]*p1; + } + } + axialVecs /= mag(axialVecs); + + + // Now we have the lambdas, curvePoints and axialVecs. + + + + info.setSize(lambdas.size()); + info = pointIndexHit(); + + // Given the current lambdas interpolate radial direction inbetween + // endpoints (all projected onto the starting coordinate system) + quaternion qStart; + vector radialStart; + { + radialStart = start-curvePoints[0]; + radialStart -= (radialStart&axialVecs[0])*axialVecs[0]; + radialStart /= mag(radialStart); + qStart = quaternion(radialStart, 0.0); + + info[0] = pointIndexHit(true, start, 0); + } + + quaternion qProjectedEnd; + { + vector radialEnd(end-curvePoints.last()); + radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last(); + radialEnd /= mag(radialEnd); + + vector projectedEnd = radialEnd; + projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0]; + projectedEnd /= mag(projectedEnd); + qProjectedEnd = quaternion(projectedEnd, 0.0); + + info.last() = pointIndexHit(true, end, 0); + } + + for (label i = 1; i < lambdas.size()-1; i++) + { + quaternion q(slerp(qStart, qProjectedEnd, lambdas[i])); + vector radialDir(q.transform(radialStart)); + + radialDir -= (radialDir&axialVecs[i])*axialVecs.last(); + radialDir /= mag(radialDir); + + info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0); + } +} + + +void Foam::searchableCurve::getRegion +( + const List<pointIndexHit>& info, + labelList& region +) const +{ + region.setSize(info.size()); + region = 0; +} + + +void Foam::searchableCurve::getNormal +( + const List<pointIndexHit>& info, + vectorField& normal +) const +{ + const edgeMesh& mesh = eMeshPtr_(); + const indexedOctree<treeDataEdge>& tree = edgeTree_(); + const edgeList& edges = mesh.edges(); + const pointField& points = mesh.points(); + + normal.setSize(info.size()); + normal = Zero; + + forAll(info, i) + { + if (info[i].hit()) + { + // Find nearest on curve + pointIndexHit curvePt = tree.findNearest + ( + info[i].hitPoint(), + Foam::magSqr(bounds().span()) + ); + + normal[i] = info[i].hitPoint()-curvePt.hitPoint(); + + // Subtract axial direction + vector axialVec = edges[curvePt.index()].vec(points); + axialVec /= mag(axialVec); + normal[i] -= (normal[i]&axialVec)*axialVec; + normal[i] /= mag(normal[i]); + } + } +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/searchableCurve/searchableCurve.H b/src/mesh/blockMesh/searchableCurve/searchableCurve.H new file mode 100644 index 0000000000000000000000000000000000000000..3c5f122cbff4653a5126ccfe4b3e1ac3b4e882b5 --- /dev/null +++ b/src/mesh/blockMesh/searchableCurve/searchableCurve.H @@ -0,0 +1,245 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::searchableCurve + +Description + Searching on edgemesh with constant radius + +SourceFiles + searchableCurve.C + +\*---------------------------------------------------------------------------*/ + +#ifndef searchableCurve_H +#define searchableCurve_H + +#include "treeBoundBox.H" +#include "searchableSurface.H" +//#include "edgeMesh.H" +//#include "indexedOctree.H" +//#include "treeDataEdge.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class edgeMesh; +class treeDataEdge; +template <class Type> class indexedOctree; + +/*---------------------------------------------------------------------------*\ + Class searchableCurve Declaration +\*---------------------------------------------------------------------------*/ + +class searchableCurve +: + public searchableSurface +{ +private: + + // Private Member Data + + //- Feature + autoPtr<edgeMesh> eMeshPtr_; + + //- Search structure + autoPtr<indexedOctree<treeDataEdge>> edgeTree_; + + //- Radius + const scalar radius_; + + //- Names of regions + mutable wordList regions_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + searchableCurve(const searchableCurve&); + + //- Disallow default bitwise assignment + void operator=(const searchableCurve&); + + +public: + + //- Runtime type information + TypeName("searchableCurve"); + + + // Constructors + + //- Construct from dictionary (used by searchableSurface) + searchableCurve + ( + const IOobject& io, + const dictionary& dict + ); + + + //- Destructor + virtual ~searchableCurve(); + + + // Member Functions + + virtual const wordList& regions() const; + + //- Whether supports volume type below + virtual bool hasVolumeType() const + { + return true; + } + + //- Range of local indices that can be returned. + virtual label size() const; + + //- Get representative set of element coordinates + // Usually the element centres (should be of length size()). + virtual tmp<pointField> coordinates() const; + + //- Get bounding spheres (centre and radius squared), one per element. + // Any point on element is guaranteed to be inside. + virtual void boundingSpheres + ( + pointField& centres, + scalarField& radiusSqr + ) const; + + //- Get the points that define the surface. + virtual tmp<pointField> points() const + { + return coordinates(); + } + + //- Does any part of the surface overlap the supplied bound box? + virtual bool overlaps(const boundBox& bb) const + { + NotImplemented; + return false; + } + + + // Multiple point queries. + + virtual void findNearest + ( + const pointField& sample, + const scalarField& nearestDistSqr, + List<pointIndexHit>& + ) const; + + //- Unique to parametric geometry: given points find + // an interpolated (along the curve) point on the surface. + // The lambdas[0] is equivalent for start, lambdas.last() + // is equivalent for end. + virtual void findNearest + ( + const point& start, + const point& end, + const scalarField& lambdas, + const scalarField& nearestDistSqr, + List<pointIndexHit>& + ) const; + + virtual void findLine + ( + const pointField& start, + const pointField& end, + List<pointIndexHit>& + ) const + { + NotImplemented; + } + + virtual void findLineAny + ( + const pointField& start, + const pointField& end, + List<pointIndexHit>& + ) const + { + NotImplemented; + } + + //- Get all intersections in order from start to end. + virtual void findLineAll + ( + const pointField& start, + const pointField& end, + List<List<pointIndexHit>>& + ) const + { + NotImplemented; + } + + //- From a set of points and indices get the region + virtual void getRegion + ( + const List<pointIndexHit>&, + labelList& region + ) const; + + //- From a set of points and indices get the normal + virtual void getNormal + ( + const List<pointIndexHit>&, + vectorField& normal + ) const; + + //- Determine type (inside/outside/mixed) for point. unknown if + // cannot be determined (e.g. non-manifold surface) + virtual void getVolumeType + ( + const pointField&, + List<volumeType>& + ) const + { + NotImplemented; + } + + + // regIOobject implementation + + bool writeData(Ostream&) const + { + NotImplemented; + return false; + } + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/pipe/Allrun b/tutorials/mesh/blockMesh/pipe/Allrun new file mode 100755 index 0000000000000000000000000000000000000000..3c4ef3e7fffd9530e7208b7e0b463a777a3bac00 --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/Allrun @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +runApplication blockMesh + +#------------------------------------------------------------------------------ diff --git a/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj new file mode 100644 index 0000000000000000000000000000000000000000..9d5e40f6d7725c20f2a5a00c5e608bbaaf46e24f --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj @@ -0,0 +1,9 @@ +v -5 0.4 0 +v -4 0.4 0 +v -3 0.4 0.7 +v -2.5 0.4 2 +v -2 0.4 2 +v -1 0.4 0.7 +v -0.6 0.4 0 +v 10 0.4 0 +l 1 2 3 4 5 6 7 8 diff --git a/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk new file mode 100644 index 0000000000000000000000000000000000000000..02142a3426c73a8d3b52d201f2221887320bb8dc --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk @@ -0,0 +1,13 @@ +# vtk DataFile Version 4.0 +vtk output +ASCII +DATASET POLYDATA +POINTS 17 double +-5.0437 0.4 -0.0159845 -4.54848 0.4 0.00946291 -4.30524 0.4 0.0219617 +-4.11475 0.4 0.0317502 -3.65463 0.4 0.0664504 -3.42509 0.4 0.242198 +-3.26981 0.4 0.570689 -3.04354 0.4 0.986036 -2.80622 0.4 1.28924 +-2.45212 0.4 1.43367 -2.10187 0.4 1.42911 -1.8115 0.4 1.2018 +-1.52708 0.4 0.866397 -1.30229 0.4 0.49514 -1.04633 0.4 0.189424 +-0.5819 0.4 -5.75752e-05 4 0.4 0 +LINES 1 18 +17 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 diff --git a/tutorials/mesh/blockMesh/pipe/system/blockMeshDict b/tutorials/mesh/blockMesh/pipe/system/blockMeshDict new file mode 100644 index 0000000000000000000000000000000000000000..39f45b9230c94e6ae647788e483013c23d012a4a --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/system/blockMeshDict @@ -0,0 +1,177 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +geometry +{ + cylinder + { + type searchableCylinder; + point1 (0 -4 0); + point2 (0 4 0); + radius 0.7; + } + cylinder3 + { + type searchableCylinder; + point1 (-10 0.4 0); + point2 (10 0.4 0); + radius 0.5; + } + cylinder2 + { + type searchableCurve; + file "curve2.vtk"; + radius 0.5; + } + inletPlane + { + type searchablePlate; + origin (-4 -50 -50); + span (0 100 100); + } +} + +vertices +( + // Vertical cylinder + named v0 project (-1 -0.1 -1) (cylinder cylinder2) + named v1 project ( 1 -0.1 -1) (cylinder) + named v2 project ( 1 0.9 -1) (cylinder) + named v3 project (-1 0.9 -1) (cylinder cylinder2) + named v4 project (-1 -0.1 1) (cylinder cylinder2) + named v5 project ( 1 -0.1 1) (cylinder) + named v6 project ( 1 0.9 1) (cylinder) + named v7 project (-1 0.9 1) (cylinder cylinder2) + + + // Horizontal cylinder + named v8 project (-4 0 -0.5) (cylinder2 inletPlane) + named v9 project (-4 1 -0.5) (cylinder2 inletPlane) + named v10 project (-4 0 0.5) (cylinder2 inletPlane) + named v11 project (-4 1 0.5) (cylinder2 inletPlane) + + + // On top of vertical cylinder + named v12 project (-1 2 -1) (cylinder) + named v13 project ( 1 2 -1) (cylinder) + named v14 project ( 1 2 1) (cylinder) + named v15 project (-1 2 1) (cylinder) + + // Below vertical cylinder + named v16 project (-1 -1 -1) (cylinder) + named v17 project ( 1 -1 -1) (cylinder) + named v18 project ( 1 -1 1) (cylinder) + named v19 project (-1 -1 1) (cylinder) +); + +blocks +( + hex (v0 v1 v2 v3 v4 v5 v6 v7) (8 8 8) simpleGrading (1 1 1) + named sideBlock hex (v0 v3 v9 v8 v4 v7 v11 v10) (8 20 8) + simpleGrading (1 1 1) + + hex ( v7 v6 v2 v3 v15 v14 v13 v12) (8 8 8) simpleGrading (1 1 1) + hex (v16 v19 v18 v17 v0 v4 v5 v1) (8 8 8) simpleGrading (1 1 1) +); + +edges +( + project v0 v1 (cylinder) + project v1 v2 (cylinder) + project v2 v3 (cylinder) + + project v1 v5 (cylinder) + project v2 v6 (cylinder) + + project v4 v5 (cylinder) + project v5 v6 (cylinder) + project v6 v7 (cylinder) + + // Common face + project v3 v0 (cylinder cylinder2) + project v3 v7 (cylinder cylinder2) + project v7 v4 (cylinder cylinder2) + project v0 v4 (cylinder cylinder2) + + // Inlet + project v8 v10 (cylinder2 inletPlane) + project v10 v11 (cylinder2 inletPlane) + project v11 v9 (cylinder2 inletPlane) + project v9 v8 (cylinder2 inletPlane) + + // Sides of horizontal cylinder. Use projectCurve to do interpolation + // for radial direction to keep points along edges at constant radial + // direction. + projectCurve v8 v0 (cylinder2) + projectCurve v9 v3 (cylinder2) + projectCurve v11 v7 (cylinder2) + projectCurve v10 v4 (cylinder2) + + + + // Top cylinder + project v12 v15 (cylinder) + project v15 v14 (cylinder) + project v14 v13 (cylinder) + project v13 v12 (cylinder) + + // Bottom cylinder + project v16 v17 (cylinder) + project v17 v18 (cylinder) + project v18 v19 (cylinder) + project v19 v16 (cylinder) +); + +faces +( + // Common face + project (v0 v4 v7 v3) cylinder + + project (v8 v0 v4 v10) cylinder2 + project (v10 v4 v7 v11) cylinder2 + project (v11 v7 v3 v9) cylinder2 + project (v8 v9 v3 v0) cylinder2 +); + + +defaultPatch +{ + name walls; + type wall; +} + +boundary +( + side + { + type patch; + faces ((sideBlock 3)); //((v8 v10 v11 v9)); + } + + inlet + { + type patch; + faces ((v17 v18 v19 v16)); + } + + outlet + { + type patch; + faces ((v12 v15 v14 v13)); + } +); + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/pipe/system/controlDict b/tutorials/mesh/blockMesh/pipe/system/controlDict new file mode 100644 index 0000000000000000000000000000000000000000..728c8b0b1c27b65296adf5b2f5a1bf7f5a7400dc --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/system/controlDict @@ -0,0 +1,58 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +libs ("libblockMesh.so"); + +DebugSwitches +{ +// project 1; +// searchableCurve 1; +// projectCurve 1; +} + +application blockMesh; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 0; + +deltaT 0; + +writeControl timeStep; + +writeInterval 1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/pipe/system/fvSchemes b/tutorials/mesh/blockMesh/pipe/system/fvSchemes new file mode 100644 index 0000000000000000000000000000000000000000..606b28ce04cbf3cc6330261345fabaa76b4d47b4 --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/system/fvSchemes @@ -0,0 +1,37 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{} + +gradSchemes +{} + +divSchemes +{} + +laplacianSchemes +{} + +interpolationSchemes +{} + +snGradSchemes +{} + + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/pipe/system/fvSolution b/tutorials/mesh/blockMesh/pipe/system/fvSolution new file mode 100644 index 0000000000000000000000000000000000000000..6db4e1a6ba541b790f69ca0a506d8c2f547e51cc --- /dev/null +++ b/tutorials/mesh/blockMesh/pipe/system/fvSolution @@ -0,0 +1,18 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// ************************************************************************* //