ENH: patchCloudSet: new sampledSet
[OpenFOAM-1.7.x.git] / src / sampling / sampledSurface / distanceSurface / distanceSurface.C
blob17ad23546e2e6a31a16df88d0db3d4b5f7990c2a
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
11     OpenFOAM is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by
13     the Free Software Foundation, either version 3 of the License, or
14     (at your option) any later version.
16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19     for more details.
21     You should have received a copy of the GNU General Public License
22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "distanceSurface.H"
27 #include "dictionary.H"
28 #include "volFields.H"
29 #include "volPointInterpolation.H"
30 #include "addToRunTimeSelectionTable.H"
31 #include "fvMesh.H"
32 #include "isoSurface.H"
33 // #include "isoSurfaceCell.H"
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 namespace Foam
39     defineTypeNameAndDebug(distanceSurface, 0);
40     addToRunTimeSelectionTable(sampledSurface, distanceSurface, word);
43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
45 void Foam::distanceSurface::createGeometry()
47     if (debug)
48     {
49         Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
50     }
52     // Clear any stored topologies
53     facesPtr_.clear();
55     // Clear derived data
56     clearGeom();
58     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
60     // Distance to cell centres
61     // ~~~~~~~~~~~~~~~~~~~~~~~~
63     cellDistancePtr_.reset
64     (
65         new volScalarField
66         (
67             IOobject
68             (
69                 "cellDistance",
70                 fvm.time().timeName(),
71                 fvm.time(),
72                 IOobject::NO_READ,
73                 IOobject::NO_WRITE,
74                 false
75             ),
76             fvm,
77             dimensionedScalar("zero", dimLength, 0)
78         )
79     );
80     volScalarField& cellDistance = cellDistancePtr_();
82     // Internal field
83     {
84         const pointField& cc = fvm.C();
85         scalarField& fld = cellDistance.internalField();
87         List<pointIndexHit> nearest;
88         surfPtr_().findNearest
89         (
90             cc,
91             scalarField(cc.size(), GREAT),
92             nearest
93         );
95         if (signed_)
96         {
97             vectorField normal;
98             surfPtr_().getNormal(nearest, normal);
100             forAll(nearest, i)
101             {
102                 vector d(cc[i]-nearest[i].hitPoint());
104                 if ((d&normal[i]) > 0)
105                 {
106                     fld[i] = Foam::mag(d);
107                 }
108                 else
109                 {
110                     fld[i] = -Foam::mag(d);
111                 }
112             }
113         }
114         else
115         {
116             forAll(nearest, i)
117             {
118                 fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
119             }
120         }
121     }
123     // Patch fields
124     {
125         forAll(fvm.C().boundaryField(), patchI)
126         {
127             const pointField& cc = fvm.C().boundaryField()[patchI];
128             fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
130             List<pointIndexHit> nearest;
131             surfPtr_().findNearest
132             (
133                 cc,
134                 scalarField(cc.size(), GREAT),
135                 nearest
136             );
138             if (signed_)
139             {
140                 vectorField normal;
141                 surfPtr_().getNormal(nearest, normal);
143                 forAll(nearest, i)
144                 {
145                     vector d(cc[i]-nearest[i].hitPoint());
147                     if ((d&normal[i]) > 0)
148                     {
149                         fld[i] = Foam::mag(d);
150                     }
151                     else
152                     {
153                         fld[i] = -Foam::mag(d);
154                     }
155                 }
156             }
157             else
158             {
159                 forAll(nearest, i)
160                 {
161                     fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
162                 }
163             }
164         }
165     }
168     // On processor patches the mesh.C() will already be the cell centre
169     // on the opposite side so no need to swap cellDistance.
172     // Distance to points
173     pointDistance_.setSize(fvm.nPoints());
174     {
175         const pointField& pts = fvm.points();
177         List<pointIndexHit> nearest;
178         surfPtr_().findNearest
179         (
180             pts,
181             scalarField(pts.size(), GREAT),
182             nearest
183         );
185         if (signed_)
186         {
187             vectorField normal;
188             surfPtr_().getNormal(nearest, normal);
190             forAll(nearest, i)
191             {
192                 vector d(pts[i]-nearest[i].hitPoint());
194                 if ((d&normal[i]) > 0)
195                 {
196                     pointDistance_[i] = Foam::mag(d);
197                 }
198                 else
199                 {
200                     pointDistance_[i] = -Foam::mag(d);
201                 }
202             }
203         }
204         else
205         {
206             forAll(nearest, i)
207             {
208                 pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
209             }
210         }
211     }
214     if (debug)
215     {
216         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
217         cellDistance.write();
218         pointScalarField pDist
219         (
220             IOobject
221             (
222                 "pointDistance",
223                 fvm.time().timeName(),
224                 fvm.time(),
225                 IOobject::NO_READ,
226                 IOobject::NO_WRITE,
227                 false
228             ),
229             pointMesh::New(fvm),
230             dimensionedScalar("zero", dimLength, 0)
231         );
232         pDist.internalField() = pointDistance_;
234         Pout<< "Writing point distance:" << pDist.objectPath() << endl;
235         pDist.write();
236     }
239     //- Direct from cell field and point field.
240     isoSurfPtr_.reset
241     (
242         new isoSurface
243         (
244             cellDistance,
245             pointDistance_,
246             distance_,
247             regularise_
248         )
249     );
251     if (debug)
252     {
253         print(Pout);
254         Pout<< endl;
255     }
259 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
261 Foam::distanceSurface::distanceSurface
263     const word& name,
264     const polyMesh& mesh,
265     const dictionary& dict
268     sampledSurface(name, mesh, dict),
269     surfPtr_
270     (
271         searchableSurface::New
272         (
273             dict.lookup("surfaceType"),
274             IOobject
275             (
276                 dict.lookupOrDefault("surfaceName", name),  // name
277                 mesh.time().constant(),                     // directory
278                 "triSurface",                               // instance
279                 mesh.time(),                                // registry
280                 IOobject::MUST_READ,
281                 IOobject::NO_WRITE
282             ),
283             dict
284         )
285     ),
286     distance_(readScalar(dict.lookup("distance"))),
287     signed_(readBool(dict.lookup("signed"))),
288     regularise_(dict.lookupOrDefault("regularise", true)),
289     zoneName_(word::null),
290     needsUpdate_(true),
291     isoSurfPtr_(NULL),
292     facesPtr_(NULL)
294 //    dict.readIfPresent("zone", zoneName_);
296 //    if (debug && zoneName_.size())
297 //    {
298 //        if (mesh.cellZones().findZoneID(zoneName_) < 0)
299 //        {
300 //            Info<< "cellZone \"" << zoneName_
301 //                << "\" not found - using entire mesh" << endl;
302 //        }
303 //    }
307 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
309 Foam::distanceSurface::~distanceSurface()
313 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
315 bool Foam::distanceSurface::needsUpdate() const
317     return needsUpdate_;
321 bool Foam::distanceSurface::expire()
323     if (debug)
324     {
325         Pout<< "distanceSurface::expire :"
326             << " have-facesPtr_:" << facesPtr_.valid()
327             << " needsUpdate_:" << needsUpdate_ << endl;
328     }
330     // Clear any stored topologies
331     facesPtr_.clear();
333     // Clear derived data
334     clearGeom();
336     // already marked as expired
337     if (needsUpdate_)
338     {
339         return false;
340     }
342     needsUpdate_ = true;
343     return true;
347 bool Foam::distanceSurface::update()
349     if (debug)
350     {
351         Pout<< "distanceSurface::update :"
352             << " have-facesPtr_:" << facesPtr_.valid()
353             << " needsUpdate_:" << needsUpdate_ << endl;
354     }
356     if (!needsUpdate_)
357     {
358         return false;
359     }
361     createGeometry();
363     needsUpdate_ = false;
364     return true;
368 Foam::tmp<Foam::scalarField>
369 Foam::distanceSurface::sample
371     const volScalarField& vField
372 ) const
374     return sampleField(vField);
378 Foam::tmp<Foam::vectorField>
379 Foam::distanceSurface::sample
381     const volVectorField& vField
382 ) const
384     return sampleField(vField);
388 Foam::tmp<Foam::sphericalTensorField>
389 Foam::distanceSurface::sample
391     const volSphericalTensorField& vField
392 ) const
394     return sampleField(vField);
398 Foam::tmp<Foam::symmTensorField>
399 Foam::distanceSurface::sample
401     const volSymmTensorField& vField
402 ) const
404     return sampleField(vField);
408 Foam::tmp<Foam::tensorField>
409 Foam::distanceSurface::sample
411     const volTensorField& vField
412 ) const
414     return sampleField(vField);
418 Foam::tmp<Foam::scalarField>
419 Foam::distanceSurface::interpolate
421     const interpolation<scalar>& interpolator
422 ) const
424     return interpolateField(interpolator);
428 Foam::tmp<Foam::vectorField>
429 Foam::distanceSurface::interpolate
431     const interpolation<vector>& interpolator
432 ) const
434     return interpolateField(interpolator);
437 Foam::tmp<Foam::sphericalTensorField>
438 Foam::distanceSurface::interpolate
440     const interpolation<sphericalTensor>& interpolator
441 ) const
443     return interpolateField(interpolator);
447 Foam::tmp<Foam::symmTensorField>
448 Foam::distanceSurface::interpolate
450     const interpolation<symmTensor>& interpolator
451 ) const
453     return interpolateField(interpolator);
457 Foam::tmp<Foam::tensorField>
458 Foam::distanceSurface::interpolate
460     const interpolation<tensor>& interpolator
461 ) const
463     return interpolateField(interpolator);
467 void Foam::distanceSurface::print(Ostream& os) const
469     os  << "distanceSurface: " << name() << " :"
470         << "  surface:" << surfPtr_().name()
471         << "  distance:" << distance_
472         << "  faces:" << faces().size()
473         << "  points:" << points().size();
477 // ************************************************************************* //