src/finiteArea/faMesh/faMeshUpdate.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "faMesh.H"
  29 #include "mapPolyMesh.H"
  30 #include "MapFaFields.H"
  31 #include "faMeshMapper.H"
  32 #include "areaFields.H"
  33 #include "edgeFields.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 bool Foam::faMesh::updateMesh(const mapPolyMesh& mpm) const
  38 {
  39     if (debug)
  40     {
  41         Info<< "bool faMesh::updateMesh(const mapPolyMesh& mpm) const : "
  42             << "Updating mesh" << endl;
  43     }
  44
  45     if (!mpm.morphing())
  46     {
  47         // No topo change
  48         return false;
  49     }
  50
  51     // Create fa mesh mapper, using the old mesh
  52     const faMeshMapper mapper(*this, mpm);
  53
  54
  55     // Rebuild mesh
  56
  57     // Cast away const for interface reasons.  HJ, 12/Aug/2011
  58     faMesh& m = const_cast<faMesh&>(*this);
  59
  60
  61     // Clear existing mesh data
  62     clearOut();
  63
  64     // Set new labels
  65     m.faceLabels_ = mapper.areaMap().newFaceLabels();
  66
  67     const indirectPrimitivePatch& bp = patch();
  68
  69     // Collect patch data
  70     const label nTotalEdges = bp.nEdges();
  71     const label nInternalEdges = bp.nInternalEdges();
  72     const labelListList& edgeFaces = bp.edgeFaces();
  73
  74     labelListList patchEdges(boundary_.size());
  75
  76     // Special handling required for faces that have more than one edge
  77     // Each patch will be visited separately
  78
  79     labelList edgeToPatch(nTotalEdges - nInternalEdges, -1);
  80     const labelList& newFaceLabelsMap = mapper.areaMap().newFaceLabelsMap();
  81
  82     const labelListList& oldPatchEdgeFaces = mapper.oldPatchEdgeFaces();
  83
  84     forAll (oldPatchEdgeFaces, patchI)
  85     {
  86         labelList& curPatchEdges = patchEdges[patchI];
  87         curPatchEdges.setSize(nTotalEdges - nInternalEdges);
  88         label nCurPatchEdges = 0;
  89
  90         // Note: it is possible to pick up the old-to-new boundary patch
  91         // mapping, but currently this is not done.  HJ, 13/Aug/2011
  92
  93         // Make a fast lookup
  94         labelHashSet oldFaceLookup(oldPatchEdgeFaces[patchI]);
  95
  96         for  (label edgeI = nInternalEdges; edgeI < nTotalEdges; edgeI++)
  97         {
  98             if (edgeToPatch[edgeI - nInternalEdges] > -1)
  99             {
 100                 // Edge already found; continue with the next one
 101                 continue;
 102             }
 103
 104             // Boundary edges will only have one face next to them
 105             const label oldFaceIndex = newFaceLabelsMap[edgeFaces[edgeI][0]];
 106
 107             if (oldFaceIndex > -1)
 108             {
 109                 // Old face exists.  See if it has got an edge in this patch
 110                 if (oldFaceLookup.found(oldFaceIndex))
 111                 {
 112                     // Face found, add it to the patch
 113                     curPatchEdges[nCurPatchEdges] = edgeI;
 114                     nCurPatchEdges++;
 115
 116                     edgeToPatch[edgeI - nInternalEdges] = patchI;
 117                 }
 118             }
 119         }
 120
 121         // Collected all faces for the current patch
 122         curPatchEdges.setSize(nCurPatchEdges);
 123     }
 124
 125     // Set new edges for all patches
 126     forAll (m.boundary_, patchI)
 127     {
 128         m.boundary_[patchI].resetEdges(patchEdges[patchI]);
 129     }
 130
 131     m.setPrimitiveMeshData();
 132
 133     // Create global mesh data
 134     if (Pstream::parRun())
 135     {
 136         globalData();
 137     }
 138
 139     // Calculate topology for the patches (processor-processor comms etc.)
 140     m.boundary_.updateMesh();
 141
 142     // Calculate the geometry for the patches (transformation tensors etc.)
 143     m.boundary_.calcGeometry();
 144
 145
 146     // Map fields
 147     mapFields(mapper);
 148
 149     // Map old areas
 150     mapOldAreas(mapper);
 151
 152     // Update edge interpolation
 153     edgeInterpolation::movePoints();
 154
 155     return true;
 156 }
 157
 158
 159 void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
 160 {
 161     // Map all the areaFields in the objectRegistry
 162     MapGeometricFields<scalar, faPatchField, faMeshMapper, areaMesh>(mapper);
 163     MapGeometricFields<vector, faPatchField, faMeshMapper, areaMesh>(mapper);
 164     MapGeometricFields<sphericalTensor, faPatchField, faMeshMapper, areaMesh>
 165         (mapper);
 166     MapGeometricFields<symmTensor, faPatchField, faMeshMapper, areaMesh>
 167         (mapper);
 168     MapGeometricFields<tensor, faPatchField, faMeshMapper, areaMesh>(mapper);
 169
 170     // Map all the edgeFields in the objectRegistry
 171     MapGeometricFields<scalar, faePatchField, faMeshMapper, edgeMesh>(mapper);
 172     MapGeometricFields<vector, faePatchField, faMeshMapper, edgeMesh>(mapper);
 173     MapGeometricFields<sphericalTensor, faePatchField, faMeshMapper, edgeMesh>
 174         (mapper);
 175     MapGeometricFields<symmTensor, faePatchField, faMeshMapper, edgeMesh>
 176         (mapper);
 177     MapGeometricFields<tensor, faePatchField, faMeshMapper, edgeMesh>(mapper);
 178 }
 179
 180
 181 void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
 182 {
 183     if (S0Ptr_)
 184     {
 185         if (debug)
 186         {
 187             InfoIn("void faMesh::mapOldAreas(const faMeshMapper& mapper)")
 188                 << "Mapping old face areas." << endl;
 189         }
 190
 191         scalarField& S0 = *S0Ptr_;
 192
 193         scalarField savedS0(S0);
 194         S0.setSize(nFaces());
 195
 196         const labelList& faceMap = mapper.areaMap().newFaceLabelsMap();
 197
 198         // Map existing old areas; for new faces set area to zero
 199         forAll (faceMap, faceI)
 200         {
 201             if (faceMap[faceI] > -1)
 202             {
 203                 S0[faceI] = savedS0[faceMap[faceI]];
 204             }
 205             else
 206             {
 207                 S0[faceI] = 0;
 208             }
 209         }
 210     }
 211
 212     if (S00Ptr_)
 213     {
 214         if (debug)
 215         {
 216             InfoIn("void faMesh::mapOldAreas(const faMeshMapper& mapper)")
 217                 << "Mapping old-old face areas." << endl;
 218         }
 219
 220         scalarField& S00 = *S00Ptr_;
 221
 222         scalarField savedS00(S00);
 223         S00.setSize(nFaces());
 224
 225         const labelList& faceMap = mapper.areaMap().newFaceLabelsMap();
 226
 227         // Map old areas for existing faces; for new faces, set area to zero
 228         forAll (faceMap, faceI)
 229         {
 230             if (faceMap[faceI] > -1)
 231             {
 232                 S00[faceI] = savedS00[faceMap[faceI]];
 233             }
 234             else
 235             {
 236                 S00[faceI] = 0;
 237             }
 238         }
 239     }
 240
 241 }
 242
 243
 244 // ************************************************************************* //