src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshSimplex.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 "Map.H"
  29 #include "partTriMeshSimplex.H"
  30
  31 //#define DEBUGSmooth
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 partTriMeshSimplex::partTriMeshSimplex
  41 (
  42     const partTriMesh& tm,
  43     const label pI
  44 )
  45 :
  46     pts_(),
  47     trias_()
  48 {
  49     const pointField& points = tm.points();
  50     const LongList<labelledTri>& trias = tm.triangles();
  51     const VRWGraph& pt = tm.pointTriangles();
  52     const LongList<direction>& pType = tm.pointType();
  53
  54     //trias_.setSize(pt.sizeOfRow(pI));
  55     label counter(0);
  56
  57     Map<label> addr(2*pt.sizeOfRow(pI));
  58     forAllRow(pt, pI, tI)
  59     {
  60         const labelledTri& tri = trias[pt(pI, tI)];
  61         for(label i=0;i<3;++i)
  62         {
  63             const label tpI = tri[i];
  64
  65             if( !addr.found(tpI) )
  66             {
  67                 addr.insert(tpI, counter);
  68                 pts_.append(points[tpI]);
  69                 ++counter;
  70             }
  71         }
  72
  73         # ifdef DEBUGSmooth
  74         Info << "Tet " << tetI << " is " << tet << endl;
  75         # endif
  76
  77         label pos(-1);
  78         for(label i=0;i<3;++i)
  79             if( tri[i] == pI )
  80             {
  81                 pos = i;
  82                 break;
  83             }
  84
  85         //- avoid using triangle serving as barriers for other points
  86         if( !(pType[tri[2]] & partTriMesh::FACECENTRE) && (pos != 0) )
  87             continue;
  88
  89         switch( pos )
  90         {
  91             case 0:
  92             {
  93                 trias_.append
  94                 (
  95                     triFace(addr[tri[0]], addr[tri[1]], addr[tri[2]])
  96                 );
  97             } break;
  98             case 1:
  99             {
 100                 trias_.append
 101                 (
 102                     triFace(addr[tri[1]], addr[tri[2]], addr[tri[0]])
 103                 );
 104             } break;
 105             case 2:
 106             {
 107                 trias_.append
 108                 (
 109                     triFace(addr[tri[2]], addr[tri[0]], addr[tri[1]])
 110                 );
 111             } break;
 112             default:
 113             {
 114                 FatalErrorIn
 115                 (
 116                     "partTriMeshSimplex::partTriMeshSimplex("
 117                     "(const partTriMesh& tm, const label pI)"
 118                 ) << "Point " << pI << " is not present in triangle" << tri
 119                     << abort(FatalError);
 120             }
 121         }
 122     }
 123
 124     # ifdef DEBUGSmooth
 125     Info << "Simplex at point " << pI << " points " << pts_ << endl;
 126     Info << "Simplex at point " << pI << " triangles " << trias_ << endl;
 127     # endif
 128
 129     if( pts_.size() == 0 || trias_.size() == 0 )
 130         FatalError << "Simplex at point " << pI << " is not valid "
 131                    << pts_ << " triangles " << trias_ << abort(FatalError);
 132 }
 133
 134 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 135
 136 partTriMeshSimplex::~partTriMeshSimplex()
 137 {}
 138
 139 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 140
 141 vector partTriMeshSimplex::normal() const
 142 {
 143     vector normal(vector::zero);
 144     scalar magN(0.0);
 145
 146     forAll(trias_, tI)
 147     {
 148         const triFace& t = trias_[tI];
 149
 150         vector n
 151         (
 152             0.5 * ((pts_[t[1]] - pts_[t[0]]) ^ (pts_[t[2]] - pts_[t[0]]))
 153         );
 154         const scalar magn = mag(n);
 155
 156         normal += n;
 157         magN += magn;
 158     }
 159
 160     return (normal / (magN + VSMALL));
 161 }
 162
 163 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 164
 165 } // End namespace Foam
 166
 167 // ************************************************************************* //