FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

primitiveMeshEdges.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "primitiveMesh.H"
00027 #include <OpenFOAM/DynamicList.H>
00028 #include <OpenFOAM/demandDrivenData.H>
00029 #include <OpenFOAM/SortableList.H>
00030 #include <OpenFOAM/ListOps.H>
00031 
00032 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00033 
00034 // Returns edgeI between two points.
00035 Foam::label Foam::primitiveMesh::getEdge
00036 (
00037     List<DynamicList<label> >& pe,
00038     DynamicList<edge>& es,
00039 
00040     const label pointI,
00041     const label nextPointI
00042 )
00043 {
00044     // Find connection between pointI and nextPointI
00045     forAll(pe[pointI], ppI)
00046     {
00047         label eI = pe[pointI][ppI];
00048 
00049         const edge& e = es[eI];
00050 
00051         if (e.start() == nextPointI || e.end() == nextPointI)
00052         {
00053             return eI;
00054         }
00055     }
00056 
00057     // Make new edge.
00058     label edgeI = es.size();
00059     pe[pointI].append(edgeI);
00060     pe[nextPointI].append(edgeI);
00061     if (pointI < nextPointI)
00062     {
00063         es.append(edge(pointI, nextPointI));
00064     }
00065     else
00066     {
00067         es.append(edge(nextPointI, pointI));
00068     }
00069     return edgeI;
00070 }
00071 
00072 
00073 void Foam::primitiveMesh::calcEdges(const bool doFaceEdges) const
00074 {
00075     if (debug)
00076     {
00077         Pout<< "primitiveMesh::calcEdges(const bool) : "
00078             << "calculating edges, pointEdges and optionally faceEdges"
00079             << endl;
00080     }
00081 
00082     // It is an error to attempt to recalculate edges
00083     // if the pointer is already set
00084     if ((edgesPtr_ || pePtr_) || (doFaceEdges && fePtr_))
00085     {
00086         FatalErrorIn("primitiveMesh::calcEdges(const bool) const")
00087             << "edges or pointEdges or faceEdges already calculated"
00088             << abort(FatalError);
00089     }
00090     else
00091     {
00092         // ALGORITHM:
00093         // Go through the faces list. Search pointEdges for existing edge.
00094         // If not found create edge and add to pointEdges.
00095         // In second loop sort edges in order of increasing neighbouring
00096         // point.
00097         // This algorithm replaces the one using pointFaces which used more
00098         // allocations but less memory and was on practical cases
00099         // quite a bit slower.
00100 
00101         const faceList& fcs = faces();
00102 
00103         // Size up lists
00104         // ~~~~~~~~~~~~~
00105 
00106         // Estimate pointEdges storage
00107         List<DynamicList<label> > pe(nPoints());
00108         forAll(pe, pointI)
00109         {
00110             pe[pointI].setCapacity(primitiveMesh::edgesPerPoint_);
00111         }
00112 
00113         // Estimate edges storage
00114         DynamicList<edge> es(pe.size()*primitiveMesh::edgesPerPoint_/2);
00115 
00116         // Estimate faceEdges storage
00117         if (doFaceEdges)
00118         {
00119             fePtr_ = new labelListList(fcs.size());
00120             labelListList& faceEdges = *fePtr_;
00121             forAll(fcs, faceI)
00122             {
00123                 faceEdges[faceI].setSize(fcs[faceI].size());
00124             }
00125         }
00126 
00127 
00128         // Find consecutive face points in edge list
00129         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00130 
00131         // Edges on boundary faces
00132         label nExtEdges = 0;
00133         // Edges using no boundary point
00134         nInternal0Edges_ = 0;
00135         // Edges using 1 boundary point
00136         label nInt1Edges = 0;
00137         // Edges using two boundary points but not on boundary face:
00138         // edges.size()-nExtEdges-nInternal0Edges_-nInt1Edges
00139 
00140         // Ordering is different if points are ordered.
00141         if (nInternalPoints_ == -1)
00142         {
00143             // No ordering. No distinction between types.
00144             forAll(fcs, faceI)
00145             {
00146                 const face& f = fcs[faceI];
00147 
00148                 forAll(f, fp)
00149                 {
00150                     label pointI = f[fp];
00151                     label nextPointI = f[f.fcIndex(fp)];
00152 
00153                     label edgeI = getEdge(pe, es, pointI, nextPointI);
00154 
00155                     if (doFaceEdges)
00156                     {
00157                         (*fePtr_)[faceI][fp] = edgeI;
00158                     }
00159                 }
00160             }
00161             // Assume all edges internal
00162             nExtEdges = 0;
00163             nInternal0Edges_ = es.size();
00164             nInt1Edges = 0;
00165         }
00166         else
00167         {
00168             // 1. Do external faces first. This creates external edges.
00169             for (label faceI = nInternalFaces_; faceI < fcs.size(); faceI++)
00170             {
00171                 const face& f = fcs[faceI];
00172 
00173                 forAll(f, fp)
00174                 {
00175                     label pointI = f[fp];
00176                     label nextPointI = f[f.fcIndex(fp)];
00177 
00178                     label oldNEdges = es.size();
00179                     label edgeI = getEdge(pe, es, pointI, nextPointI);
00180 
00181                     if (es.size() > oldNEdges)
00182                     {
00183                         nExtEdges++;
00184                     }
00185                     if (doFaceEdges)
00186                     {
00187                         (*fePtr_)[faceI][fp] = edgeI;
00188                     }
00189                 }
00190             }
00191 
00192             // 2. Do internal faces. This creates internal edges.
00193             for (label faceI = 0; faceI < nInternalFaces_; faceI++)
00194             {
00195                 const face& f = fcs[faceI];
00196 
00197                 forAll(f, fp)
00198                 {
00199                     label pointI = f[fp];
00200                     label nextPointI = f[f.fcIndex(fp)];
00201 
00202                     label oldNEdges = es.size();
00203                     label edgeI = getEdge(pe, es, pointI, nextPointI);
00204 
00205                     if (es.size() > oldNEdges)
00206                     {
00207                         if (pointI < nInternalPoints_)
00208                         {
00209                             if (nextPointI < nInternalPoints_)
00210                             {
00211                                 nInternal0Edges_++;
00212                             }
00213                             else
00214                             {
00215                                 nInt1Edges++;
00216                             }
00217                         }
00218                         else
00219                         {
00220                             if (nextPointI < nInternalPoints_)
00221                             {
00222                                 nInt1Edges++;
00223                             }
00224                             else
00225                             {
00226                                 // Internal edge with two points on boundary
00227                             }
00228                         }
00229                     }
00230                     if (doFaceEdges)
00231                     {
00232                         (*fePtr_)[faceI][fp] = edgeI;
00233                     }
00234                 }
00235             }
00236         }
00237 
00238 
00239         // For unsorted meshes the edges will be in order of occurrence inside
00240         // the faces. For sorted meshes the first nExtEdges will be the external
00241         // edges.
00242 
00243         if (nInternalPoints_ != -1)
00244         {
00245             nInternalEdges_ = es.size()-nExtEdges;
00246             nInternal1Edges_ = nInternal0Edges_+nInt1Edges;
00247 
00248             //Pout<< "Edge overview:" << nl
00249             //    << "    total number of edges           : " << es.size()
00250             //    << nl
00251             //    << "    boundary edges                  : " << nExtEdges
00252             //    << nl
00253             //    << "    edges using no boundary point   : "
00254             //    << nInternal0Edges_
00255             //    << nl
00256             //    << "    edges using one boundary points : " << nInt1Edges
00257             //   << nl
00258             //    << "    edges using two boundary points : "
00259             //    << es.size()-nExtEdges-nInternal0Edges_-nInt1Edges << endl;
00260         }
00261 
00262 
00263         // Like faces sort edges in order of increasing neigbouring point.
00264         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00265         // Automatically if points are sorted into internal and external points
00266         // the edges will be sorted into
00267         // - internal point to internal point
00268         // - internal point to external point
00269         // - external point to external point
00270         // The problem is that the last one mixes external edges with internal
00271         // edges with two points on the boundary.
00272 
00273 
00274         // Map to sort into new upper-triangular order
00275         labelList oldToNew(es.size());
00276         if (debug)
00277         {
00278             oldToNew = -1;
00279         }
00280 
00281         // start of edges with 0 boundary points
00282         label internal0EdgeI = 0;
00283 
00284         // start of edges with 1 boundary points
00285         label internal1EdgeI = nInternal0Edges_;
00286 
00287         // start of edges with 2 boundary points
00288         label internal2EdgeI = nInternal1Edges_;
00289 
00290         // start of external edges
00291         label externalEdgeI = nInternalEdges_;
00292 
00293 
00294         // To sort neighbouring points in increasing order. Defined outside
00295         // for optimisation reasons: if all connectivity size same will need
00296         // no reallocations
00297         SortableList<label> nbrPoints(primitiveMesh::edgesPerPoint_);
00298 
00299         forAll(pe, pointI)
00300         {
00301             const DynamicList<label>& pEdges = pe[pointI];
00302 
00303             nbrPoints.setSize(pEdges.size());
00304 
00305             forAll(pEdges, i)
00306             {
00307                 const edge& e = es[pEdges[i]];
00308 
00309                 label nbrPointI = e.otherVertex(pointI);
00310 
00311                 if (nbrPointI < pointI)
00312                 {
00313                     nbrPoints[i] = -1;
00314                 }
00315                 else
00316                 {
00317                     nbrPoints[i] = nbrPointI;
00318                 }
00319             }
00320             nbrPoints.sort();
00321 
00322 
00323             if (nInternalPoints_ == -1)
00324             {
00325                 // Sort all upper-triangular
00326                 forAll(nbrPoints, i)
00327                 {
00328                     if (nbrPoints[i] != -1)
00329                     {
00330                         label edgeI = pEdges[nbrPoints.indices()[i]];
00331 
00332                         oldToNew[edgeI] = internal0EdgeI++;
00333                     }
00334                 }
00335             }
00336             else
00337             {
00338                 if (pointI < nInternalPoints_)
00339                 {
00340                     forAll(nbrPoints, i)
00341                     {
00342                         label nbrPointI = nbrPoints[i];
00343 
00344                         label edgeI = pEdges[nbrPoints.indices()[i]];
00345 
00346                         if (nbrPointI != -1)
00347                         {
00348                             if (edgeI < nExtEdges)
00349                             {
00350                                 // External edge
00351                                 oldToNew[edgeI] = externalEdgeI++;
00352                             }
00353                             else if (nbrPointI < nInternalPoints_)
00354                             {
00355                                 // Both points inside
00356                                 oldToNew[edgeI] = internal0EdgeI++;
00357                             }
00358                             else
00359                             {
00360                                 // One points inside, one outside
00361                                 oldToNew[edgeI] = internal1EdgeI++;
00362                             }
00363                         }
00364                     }
00365                 }
00366                 else
00367                 {
00368                     forAll(nbrPoints, i)
00369                     {
00370                         label nbrPointI = nbrPoints[i];
00371 
00372                         label edgeI = pEdges[nbrPoints.indices()[i]];
00373 
00374                         if (nbrPointI != -1)
00375                         {
00376                             if (edgeI < nExtEdges)
00377                             {
00378                                 // External edge
00379                                 oldToNew[edgeI] = externalEdgeI++;
00380                             }
00381                             else if (nbrPointI < nInternalPoints_)
00382                             {
00383                                 // Not possible!
00384                                 FatalErrorIn("primitiveMesh::calcEdges(..)")
00385                                     << abort(FatalError);
00386                             }
00387                             else
00388                             {
00389                                 // Both points outside
00390                                 oldToNew[edgeI] = internal2EdgeI++;
00391                             }
00392                         }
00393                     }
00394                 }
00395             }
00396         }
00397 
00398 
00399         if (debug)
00400         {
00401             label edgeI = findIndex(oldToNew, -1);
00402 
00403             if (edgeI != -1)
00404             {
00405                 const edge& e = es[edgeI];
00406 
00407                 FatalErrorIn("primitiveMesh::calcEdges(..)")
00408                     << "Did not sort edge " << edgeI << " points:" << e
00409                     << " coords:" << points()[e[0]] << points()[e[1]]
00410                     << endl
00411                     << "Current buckets:" << endl
00412                     << "    internal0EdgeI:" << internal0EdgeI << endl
00413                     << "    internal1EdgeI:" << internal1EdgeI << endl
00414                     << "    internal2EdgeI:" << internal2EdgeI << endl
00415                     << "    externalEdgeI:" << externalEdgeI << endl
00416                     << abort(FatalError);
00417             }
00418         }
00419 
00420 
00421 
00422         // Renumber and transfer.
00423 
00424         // Edges
00425         edgesPtr_ = new edgeList(es.size());
00426         edgeList& edges = *edgesPtr_;
00427         forAll(es, edgeI)
00428         {
00429             edges[oldToNew[edgeI]] = es[edgeI];
00430         }
00431 
00432         // pointEdges
00433         pePtr_ = new labelListList(nPoints());
00434         labelListList& pointEdges = *pePtr_;
00435         forAll(pe, pointI)
00436         {
00437             DynamicList<label>& pEdges = pe[pointI];
00438             pEdges.shrink();
00439             inplaceRenumber(oldToNew, pEdges);
00440             pointEdges[pointI].transfer(pEdges);
00441             Foam::sort(pointEdges[pointI]);
00442         }
00443 
00444         // faceEdges
00445         if (doFaceEdges)
00446         {
00447             labelListList& faceEdges = *fePtr_;
00448             forAll(faceEdges, faceI)
00449             {
00450                 inplaceRenumber(oldToNew, faceEdges[faceI]);
00451             }
00452         }
00453     }
00454 }
00455 
00456 
00457 Foam::label Foam::primitiveMesh::findFirstCommonElementFromSortedLists
00458 (
00459     const labelList& list1,
00460     const labelList& list2
00461 )
00462 {
00463     label result = -1;
00464 
00465     labelList::const_iterator iter1 = list1.begin();
00466     labelList::const_iterator iter2 = list2.begin();
00467 
00468     while (iter1 != list1.end() && iter2 != list2.end())
00469     {
00470         if( *iter1 < *iter2)
00471         {
00472             ++iter1;
00473         }
00474         else if (*iter1 > *iter2)
00475         {
00476             ++iter2;
00477         }
00478         else
00479         {
00480             result = *iter1;
00481             break;
00482         }
00483     }
00484     if (result == -1)
00485     {
00486         FatalErrorIn
00487         (
00488             "primitiveMesh::findFirstCommonElementFromSortedLists"
00489             "(const labelList&, const labelList&)"
00490         )   << "No common elements in lists " << list1 << " and " << list2
00491             << abort(FatalError);
00492     }
00493     return result;
00494 }
00495 
00496 
00497 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00498 
00499 const Foam::edgeList& Foam::primitiveMesh::edges() const
00500 {
00501     if (!edgesPtr_)
00502     {
00503         //calcEdges(true);
00504         calcEdges(false);
00505     }
00506 
00507     return *edgesPtr_;
00508 }
00509 
00510 const Foam::labelListList& Foam::primitiveMesh::pointEdges() const
00511 {
00512     if (!pePtr_)
00513     {
00514         //calcEdges(true);
00515         calcEdges(false);
00516     }
00517 
00518     return *pePtr_;
00519 }
00520 
00521 
00522 const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
00523 {
00524     if (!fePtr_)
00525     {
00526         if (debug)
00527         {
00528             Pout<< "primitiveMesh::faceEdges() : "
00529                 << "calculating faceEdges" << endl;
00530         }
00531 
00532         //calcEdges(true);
00533         const faceList& fcs = faces();
00534         const labelListList& pe = pointEdges();
00535         const edgeList& es = edges();
00536 
00537         fePtr_ = new labelListList(fcs.size());
00538         labelListList& faceEdges = *fePtr_;
00539 
00540         forAll(fcs, faceI)
00541         {
00542             const face& f = fcs[faceI];
00543 
00544             labelList& fEdges = faceEdges[faceI];
00545             fEdges.setSize(f.size());
00546 
00547             forAll(f, fp)
00548             {
00549                 label pointI = f[fp];
00550                 label nextPointI = f[f.fcIndex(fp)];
00551 
00552                 // Find edge between pointI, nextPontI
00553                 const labelList& pEdges = pe[pointI];
00554 
00555                 forAll(pEdges, i)
00556                 {
00557                     label edgeI = pEdges[i];
00558 
00559                     if (es[edgeI].otherVertex(pointI) == nextPointI)
00560                     {
00561                         fEdges[fp] = edgeI;
00562                         break;
00563                     }
00564                 }
00565             }
00566         }
00567     }
00568 
00569     return *fePtr_;
00570 }
00571 
00572 
00573 void Foam::primitiveMesh::clearOutEdges()
00574 {
00575     deleteDemandDrivenData(edgesPtr_);
00576     deleteDemandDrivenData(pePtr_);
00577     deleteDemandDrivenData(fePtr_);
00578     labels_.clear();
00579     labelSet_.clear();
00580 }
00581 
00582 
00583 const Foam::labelList& Foam::primitiveMesh::faceEdges
00584 (
00585     const label faceI,
00586     DynamicList<label>& storage
00587 ) const
00588 {
00589     if (hasFaceEdges())
00590     {
00591         return faceEdges()[faceI];
00592     }
00593     else
00594     {
00595         const labelListList& pointEs = pointEdges();
00596         const face& f = faces()[faceI];
00597 
00598         storage.clear();
00599         if (f.size() > storage.capacity())
00600         {
00601             storage.setCapacity(f.size());
00602         }
00603 
00604         forAll(f, fp)
00605         {
00606             storage.append
00607             (
00608                 findFirstCommonElementFromSortedLists
00609                 (
00610                     pointEs[f[fp]],
00611                     pointEs[f.nextLabel(fp)]
00612                 )
00613             );
00614         }
00615 
00616         return storage;
00617     }
00618 }
00619 
00620 
00621 const Foam::labelList& Foam::primitiveMesh::faceEdges(const label faceI) const
00622 {
00623     return faceEdges(faceI, labels_);
00624 }
00625 
00626 
00627 const Foam::labelList& Foam::primitiveMesh::cellEdges
00628 (
00629     const label cellI,
00630     DynamicList<label>& storage
00631 ) const
00632 {
00633     if (hasCellEdges())
00634     {
00635         return cellEdges()[cellI];
00636     }
00637     else
00638     {
00639         const labelList& cFaces = cells()[cellI];
00640 
00641         labelSet_.clear();
00642 
00643         forAll(cFaces, i)
00644         {
00645             const labelList& fe = faceEdges(cFaces[i]);
00646 
00647             forAll(fe, feI)
00648             {
00649                 labelSet_.insert(fe[feI]);
00650             }
00651         }
00652 
00653         storage.clear();
00654 
00655         if (labelSet_.size() > storage.capacity())
00656         {
00657             storage.setCapacity(labelSet_.size());
00658         }
00659 
00660         forAllConstIter(labelHashSet, labelSet_, iter)
00661         {
00662             storage.append(iter.key());
00663         }
00664 
00665         return storage;
00666     }
00667 }
00668 
00669 
00670 const Foam::labelList& Foam::primitiveMesh::cellEdges(const label cellI) const
00671 {
00672     return cellEdges(cellI, labels_);
00673 }
00674 
00675 
00676 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00677 
00678 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines