foamVtkFvMeshAdaptorFieldTemplates.C 8.68 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
     \\/     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/>.

\*---------------------------------------------------------------------------*/

#ifndef foamVtkFvMeshAdaptorFieldTemplates_C
#define foamVtkFvMeshAdaptorFieldTemplates_C

// OpenFOAM includes
#include "error.H"
#include "emptyFvPatchField.H"
#include "wallPolyPatch.H"
#include "volPointInterpolation.H"
#include "zeroGradientFvPatchField.H"
#include "interpolatePointToCell.H"

// vtk includes
#include "vtkFloatArray.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//
// volume-fields
//

template<class Type>
void Foam::vtk::fvMeshAdaptor::convertVolField
(
    const PtrList<patchInterpolator>& patchInterpList,
    const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
    const fvMesh& mesh = fld.mesh();
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    // Interpolated field (demand driven)
    autoPtr<GeometricField<Type, pointPatchField, pointMesh>> ptfPtr;
60
    if (interpFields_)
61 62 63 64 65 66 67
    {
        ptfPtr.reset
        (
            volPointInterpolation::New(mesh).interpolate(fld).ptr()
        );
    }

68
    // INTERNAL
69 70
    convertVolFieldInternal(fld, ptfPtr);

71
    // BOUNDARY
72
    for (const label patchId : patchIds_)
73 74 75 76 77
    {
        const polyPatch& pp = patches[patchId];
        const word& longName = pp.name();

        auto iter = cachedVtp_.find(longName);
78
        if (!iter.found() || !iter.val().dataset)
79 80
        {
            // Should not happen, but for safety require a vtk geometry
81
            Pout<<"Cache miss for VTP patch " << longName << endl;
82 83 84
            continue;
        }

85
        foamVtpData& vtpData = iter.val();
86 87 88 89 90 91 92 93 94 95 96
        auto dataset = vtpData.dataset;

        // This is slightly roundabout, but we deal with groups and with
        // single patches.
        // For groups (spanning several patches) it is fairly messy to
        // get interpolated point fields. We would need to create a indirect
        // patch each time to obtain the mesh points. We thus skip that.
        //
        // To improve code reuse, we allocate the CellData as a zeroed-field
        // ahead of time.

97
        vtkSmartPointer<vtkFloatArray> cdata = vtk::Tools::zeroField<Type>
98 99 100 101 102 103 104
        (
            fld.name(),
            dataset->GetNumberOfPolys()
        );

        vtkSmartPointer<vtkFloatArray> pdata;

105
        const fvPatchField<Type>& ptf = fld.boundaryField()[patchId];
106

107 108 109 110
        if
        (
            isType<emptyFvPatchField<Type>>(ptf)
         ||
111
            (
112
                extrapPatches_
113
             && !polyPatch::constraintType(patches[patchId].type())
114
            )
115 116 117
        )
        {
            fvPatch p(ptf.patch().patch(), mesh.boundary());
118

119 120 121 122
            tmp<Field<Type>> tpptf
            (
                fvPatchField<Type>(p, fld).patchInternalField()
            );
123

124
            vtk::Tools::transcribeFloatData(cdata, tpptf());
125

126 127 128 129 130
            if
            (
                patchId < patchInterpList.size()
             && patchInterpList.set(patchId)
            )
131 132 133 134 135 136
            {
                pdata = vtk::Tools::convertFieldToVTK
                (
                    fld.name(),
                    patchInterpList[patchId].faceToPointInterpolate(tpptf)()
                );
137
            }
138 139 140
        }
        else
        {
141
            vtk::Tools::transcribeFloatData(cdata, ptf);
142

143 144 145 146 147
            if
            (
                patchId < patchInterpList.size()
             && patchInterpList.set(patchId)
            )
148
            {
149 150 151 152 153
                pdata = vtk::Tools::convertFieldToVTK
                (
                    fld.name(),
                    patchInterpList[patchId].faceToPointInterpolate(ptf)()
                );
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
            }
        }

        if (cdata)
        {
            dataset->GetCellData()->AddArray(cdata);
        }
        if (pdata)
        {
            dataset->GetPointData()->AddArray(pdata);
        }
    }
}


template<class Type>
void Foam::vtk::fvMeshAdaptor::convertVolFields
(
    const PtrList<patchInterpolator>& patchInterpList,
    const wordRes& selectFields
)
{
176
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
177 178 179

    // Restrict to GeometricField<Type, ...>

180 181 182 183 184
    for
    (
        const word& fieldName
      : mesh_.sortedNames<fieldType>(selectFields)
    )
185 186 187 188
    {
        convertVolField
        (
            patchInterpList,
189
            mesh_.lookupObject<fieldType>(fieldName)
190 191 192 193 194 195 196 197 198 199 200 201
        );
    }
}


template<class Type>
void Foam::vtk::fvMeshAdaptor::convertVolFieldInternal
(
    const GeometricField<Type, fvPatchField, volMesh>& fld,
    autoPtr<GeometricField<Type, pointPatchField, pointMesh>>& ptfPtr
)
{
202
    if (!usingInternal())
203 204 205 206
    {
        return;
    }

207
    const auto& longName = internalName();
208 209

    auto iter = cachedVtu_.find(longName);
210
    if (!iter.found() || !iter.val().dataset)
211 212
    {
        // Should not happen, but for safety require a vtk geometry
213
        Pout<<"Cache miss for VTU " << longName << endl;
214 215
        return;
    }
216
    foamVtuData& vtuData = iter.val();
217 218
    auto dataset = vtuData.dataset;

219
    dataset->GetCellData()->AddArray
220
    (
221
        vtuData.convertField(fld)
222 223 224 225
    );

    if (ptfPtr.valid())
    {
226
        dataset->GetPointData()->AddArray
227
        (
228
            convertPointField(ptfPtr(), fld, vtuData)
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
        );
    }
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<class Type>
vtkSmartPointer<vtkFloatArray> Foam::vtk::fvMeshAdaptor::convertPointField
(
    const GeometricField<Type, pointPatchField, pointMesh>& pfld,
    const GeometricField<Type, fvPatchField, volMesh>& vfld,
    const foamVtuData& vtuData
)
{
    const int nComp(pTraits<Type>::nComponents);
    const labelUList& addPointCellLabels = vtuData.additionalIds();
    const labelUList& pointMap = vtuData.pointMap();

    // Use a pointMap or address directly into mesh
    const label nPoints = (pointMap.size() ? pointMap.size() : pfld.size());

    auto data = vtkSmartPointer<vtkFloatArray>::New();
    data->SetNumberOfComponents(nComp);
    data->SetNumberOfTuples(nPoints + addPointCellLabels.size());

    // Note: using the name of the original volField
    // not the name generated by the interpolation "volPointInterpolate(<name>)"

    if (notNull(vfld))
    {
        data->SetName(vfld.name().c_str());
    }
    else
    {
        data->SetName(pfld.name().c_str());
    }

267 268 269 270 271
    DebugInfo
        << "convert Point field: " << pfld.name()
        << " size="  << (nPoints + addPointCellLabels.size())
        << " (" << nPoints << " + " << addPointCellLabels.size()
        << ") nComp=" << nComp << endl;
272 273


274 275 276
    float scratch[pTraits<Type>::nComponents];

    vtkIdType pointi = 0;
277 278
    if (pointMap.size())
    {
279
        for (const label meshPointi : pointMap)
280
        {
281
            vtk::Tools::foamToVtkTuple(scratch, pfld[meshPointi]);
282
            data->SetTuple(pointi++, scratch);
283 284 285 286
        }
    }
    else
    {
287
        for (const Type& val : pfld)
288
        {
289
            vtk::Tools::foamToVtkTuple(scratch, val);
290
            data->SetTuple(pointi++, scratch);
291 292 293 294
        }
    }

    // Continue with additional points
295
    // - correspond to cell centres
296 297 298

    if (notNull(vfld))
    {
299
        for (const label meshCelli : addPointCellLabels)
300
        {
301
            vtk::Tools::foamToVtkTuple(scratch, vfld[meshCelli]);
302
            data->SetTuple(pointi++, scratch);
303 304 305 306
        }
    }
    else
    {
307
        for (const label meshCelli : addPointCellLabels)
308
        {
309 310 311 312
            vtk::Tools::foamToVtkTuple
            (
                scratch, interpolatePointToCell(pfld, meshCelli)
            );
313
            data->SetTuple(pointi++, scratch);
314 315 316 317 318 319 320 321 322 323
        }
    }

    return data;
}


#endif

// ************************************************************************* //