src/mesh/cfMesh/utilities/faceDecomposition/faceDecomposition.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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     cfMesh 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.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "faceDecomposition.H"
  29 #include "pointField.H"
  30 #include "boolList.H"
  31 #include "helperFunctions.H"
  32
  33 // #define DEBUG_faceDecomposition
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 label faceDecomposition::concaveVertex() const
  41 {
  42     vector n = f_.normal(points_);
  43     n /= mag(n);
  44
  45     const edgeList edges = f_.edges();
  46
  47     label concaveVrt(-1);
  48
  49     forAll(edges, eI)
  50     {
  51         vector ev = edges[eI].vec(points_);
  52         ev /= mag(ev);
  53
  54         const short next = (eI+1) % f_.size();
  55         vector evn = edges[next].vec(points_);
  56         evn /= mag(evn);
  57
  58         const vector prod = (ev ^ evn);
  59
  60         if( (prod & n) < -SMALL )
  61         {
  62             if( concaveVrt != -1 )
  63             {
  64                 FatalErrorIn
  65                 (
  66                     "label faceDecomposition::concaveVertex(const label faceI) const"
  67                 ) << "Face " << f_ << " has more than one concave vertex."
  68                     << " Cannot continue ..." << exit(FatalError);
  69             }
  70
  71             label vrtIndex = edges[eI].commonVertex(edges[next]);
  72             /*
  73             if(
  74                 pointTriIndex_[vrtIndex] &&
  75                 pointTriIndex_[vrtIndex]->size() > 1
  76             )
  77             */
  78             concaveVrt = vrtIndex;
  79         }
  80     }
  81
  82     return concaveVrt;
  83 }
  84
  85 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  86 //- Constructor
  87 faceDecomposition::faceDecomposition
  88 (
  89     const face& f,
  90     const pointField& pts
  91 )
  92 :
  93     f_(f),
  94     points_(pts)
  95 {
  96 }
  97
  98 //- Destructor
  99 faceDecomposition::~faceDecomposition()
 100 {
 101 }
 102
 103 bool faceDecomposition::isFaceConvex() const
 104 {
 105     if( concaveVertex() == -1 )
 106         return true;
 107
 108     return false;
 109 }
 110
 111 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 112
 113 bool faceDecomposition::isFacePlanar(const scalar tol) const
 114 {
 115     vector nref = f_.normal(points_);
 116     nref /= mag(nref);
 117
 118     forAll(f_, pI)
 119         if( mag((points_[f_[pI]] - points_[f_[0]]) & nref) > tol )
 120         {
 121             # ifdef DEBUG_faceDecomposition
 122             Info << "Face is not planar " << endl;
 123             # endif
 124             return false;
 125         }
 126
 127     return true;
 128 }
 129
 130 bool faceDecomposition::isFacePlanar() const
 131 {
 132     scalar tol(0.0);
 133
 134     const point c = f_.centre(points_);
 135     forAll(f_, pI)
 136         tol = Foam::max(tol, Foam::mag(c - points_[f_[pI]]));
 137
 138     tol *= 0.05;
 139
 140     return isFacePlanar(tol);
 141 }
 142
 143 faceList faceDecomposition::decomposeFace() const
 144 {
 145     faceList ff = decomposeFaceIntoTriangles();
 146
 147     if( ff.size() <= 3 )
 148     {
 149         //- face is decomposed into 2 or 3 triangles. I am happy with that.
 150         return ff;
 151     }
 152
 153     boolList mergedFaces(ff.size(), false);
 154     face lf = ff[0];
 155     face rf = ff[ff.size()-1];
 156
 157     direction il(1), ir(ff.size()-2);
 158     faceList storage(2);
 159     direction fI(0);
 160
 161     while( il < ir )
 162     {
 163         //- merge on the side of lower labels
 164         face fl = help::mergeTwoFaces(ff[il], lf);
 165         if( faceDecomposition(fl, points_).isFaceConvex() )
 166         {
 167             lf = fl;
 168             if( il == ir - 1 )
 169                 storage.newElmt(fI++) = lf;
 170         }
 171         else
 172         {
 173             storage.newElmt(fI++) = lf;
 174             lf = ff[il];
 175             if( il == ir - 1 )
 176                 storage.newElmt(fI++) = lf;
 177         }
 178
 179         //- merge from the side of higher labels
 180         face fr = help::mergeTwoFaces(rf, ff[ir]);
 181         if( faceDecomposition(fr, points_).isFaceConvex() )
 182         {
 183             rf = fr;
 184             if( il == ir - 1 )
 185                 storage.newElmt(fI++) = rf;
 186         }
 187         else
 188         {
 189             storage.newElmt(fI++) = rf;
 190             rf = ff[ir];
 191             if( il == ir - 1 )
 192                 storage.newElmt(fI++) = rf;
 193         }
 194
 195         il++;
 196         ir--;
 197
 198         //- this happens if the face has odd number of edges
 199         if( il == ir )
 200         {
 201             fl = help::mergeTwoFaces(ff[il], lf);
 202             fr = help::mergeTwoFaces(rf, ff[ir]);
 203             if( faceDecomposition(fl, points_).isFaceConvex() )
 204             {
 205                 storage.newElmt(fI++) = fl;
 206                 storage.newElmt(fI++) = rf;
 207             }
 208             else if( faceDecomposition(fr, points_).isFaceConvex() )
 209             {
 210                 storage.newElmt(fI++) = fr;
 211                 storage.newElmt(fI++) = lf;
 212             }
 213             else
 214             {
 215                 storage.newElmt(fI++) = lf;
 216                 storage.newElmt(fI++) = ff[il];
 217                 storage.newElmt(fI++) = rf;
 218             }
 219         }
 220     }
 221
 222     if( storage.size() > 2 )
 223         storage.setSize(fI);
 224
 225     # ifdef DEBUG_faceDecomposition
 226     Info << "Original face " << f_ << endl;
 227     Info << "Triangles " << ff << endl;
 228     Info << "Concave vertex " << concaveVertex(f_) << endl;
 229     Info << "Storage " << storage << endl;
 230     # endif
 231
 232     return storage;
 233 }
 234
 235 faceList faceDecomposition::decomposeFaceIntoTriangles(const label cv) const
 236 {
 237     if( cv != -1 )
 238     {
 239         short start(0);
 240         forAll(f_, pI)
 241             if( cv == f_[pI] )
 242             {
 243                 start = pI;
 244                 break;
 245             }
 246
 247         faceList fcs(10);
 248         short fI(0);
 249
 250         const edgeList edg = f_.edges();
 251
 252         for(short eI=1;eI<edg.size()-1;eI++)
 253         {
 254             const short i = (eI+start) % f_.size();
 255             face add(3);
 256             add[0] = f_[start];
 257             add[1] = edg[i].start();
 258             add[2] = edg[i].end();
 259
 260             fcs.newElmt(fI++) = add;
 261         }
 262
 263         fcs.setSize(fI);
 264
 265         # ifdef DEBUG_faceDecomposition
 266         Info << "face " << faceNo << "  " << f_
 267             << " is decomposed into " << fcs << endl;
 268         # endif
 269
 270         return fcs;
 271     }
 272
 273     faceList fcs(1, f_);
 274
 275     return fcs;
 276 }
 277
 278 faceList faceDecomposition::decomposeFaceIntoTriangles() const
 279 {
 280     const label cv = concaveVertex();
 281     return decomposeFaceIntoTriangles(cv);
 282 }
 283
 284 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 285
 286 } // End namespace Foam
 287
 288 // ************************************************************************* //