Mesh.cpp 10 KB
Newer Older
Michael Turner's avatar
Michael Turner committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
//  File: Mesh.cpp
Michael Turner's avatar
Michael Turner committed
4 5 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
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
Michael Turner's avatar
Michael Turner committed
32
//  Description: Mesh object.
Michael Turner's avatar
Michael Turner committed
33 34 35
//
////////////////////////////////////////////////////////////////////////////////

Michael Turner's avatar
Michael Turner committed
36
#include <LibUtilities/BasicUtils/Progressbar.hpp>
Michael Turner's avatar
Michael Turner committed
37 38
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <NekMeshUtils/MeshElements/Mesh.h>
39

Michael Turner's avatar
Michael Turner committed
40 41 42 43
using namespace std;

namespace Nektar
{
44
namespace NekMeshUtils
Michael Turner's avatar
Michael Turner committed
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
{

/**
 * @brief Return the number of elements of the expansion dimension.
 */
unsigned int Mesh::GetNumElements()
{
    return m_element[m_expDim].size();
}

/**
 * @brief Return the number of boundary elements (i.e. one below the
 * expansion dimension).
 */
unsigned int Mesh::GetNumBndryElements()
{
    unsigned int i, nElmt = 0;

    for (i = 0; i < m_expDim; ++i)
        nElmt += m_element[i].size();

    return nElmt;
}

/**
 * @brief Return the total number of entities in the mesh (i.e. all
 * elements, regardless of dimension).
 */
unsigned int Mesh::GetNumEntities()
{
    unsigned int nEnt = 0;

    for (unsigned int d = 0; d <= m_expDim; ++d)
    {
        nEnt += m_element[d].size();
    }

    return nEnt;
}
84 85

/**
Dave Moxey's avatar
Dave Moxey committed
86 87 88 89 90 91 92 93 94 95 96 97
 * @brief Convert this mesh into a mesh of uniform polynomial order @p order
 * with a curve point distribution @p distType.
 *
 * This routine adds curvature points into a mesh so that the resulting elements
 * are all of a uniform order @p order and all high-order vertices are
 * consistently ordered. It proceeds in a bottom-up fashion:
 *
 * - First construct all edge, face and elemental geometry mappings.
 * - Then call the local MakeOrder functions on each edge, face and element of
 *   dimension Mesh::m_expDim.
 * - Finally, any boundary elements are updated so that they have the same
 *   interior degrees of freedom as their corresponding edge or face links.
98
 */
Michael Turner's avatar
Michael Turner committed
99
void Mesh::MakeOrder(int order, LibUtilities::PointsType distType)
100
{
101 102 103 104
    // Going to make a copy of the curavture information, since this is cheaper
    // than using Nektar's Geometry objects. Currently, the geometry objects
    // which make up a 3D element dont use the volume nodes, they are just
    // stored, so we can get away without copying them.
Michael Turner's avatar
Michael Turner committed
105

106 107 108 109 110
    int id = m_vertexSet.size();

    EdgeSet::iterator eit;
    FaceSet::iterator fit;

Michael Turner's avatar
Michael Turner committed
111 112
    boost::unordered_map<int, EdgeSharedPtr> edgeCopies;
    boost::unordered_map<int, FaceSharedPtr> faceCopies;
113

114 115
    // Decide on distribution of points to use for each shape type based on the
    // input we've been supplied.
116 117 118
    std::map<LibUtilities::ShapeType, LibUtilities::PointsType> pTypes;
    if (distType == LibUtilities::ePolyEvenlySpaced)
    {
119 120
        pTypes[LibUtilities::eSegment]  = LibUtilities::ePolyEvenlySpaced;
        pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced;
Dave Moxey's avatar
Dave Moxey committed
121
        pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced;
122 123 124 125
        pTypes[LibUtilities::eTetrahedron] =
            LibUtilities::eNodalTetEvenlySpaced;
        pTypes[LibUtilities::ePrism] = LibUtilities::eNodalPrismEvenlySpaced;
        pTypes[LibUtilities::eHexahedron] = LibUtilities::ePolyEvenlySpaced;
126 127 128
    }
    else if (distType == LibUtilities::eGaussLobattoLegendre)
    {
129 130
        pTypes[LibUtilities::eSegment]  = LibUtilities::eGaussLobattoLegendre;
        pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriElec;
131 132
        pTypes[LibUtilities::eQuadrilateral] =
            LibUtilities::eGaussLobattoLegendre;
Michael Turner's avatar
Michael Turner committed
133
        pTypes[LibUtilities::ePrism]       = LibUtilities::eNodalPrismElec;
134 135
        pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec;
        pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre;
136
    }
Dave Moxey's avatar
Dave Moxey committed
137 138 139 140
    else
    {
        ASSERTL1(false, "Mesh::MakeOrder does not support this points type.");
    }
141

142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
    // Begin by copying mesh objects for edges and faces so that we don't affect
    // any neighbouring elements in the mesh as we process each element. At the
    // same time we delete the curvature from the original edge and face, which
    // will be re-added with the MakeOrder routine.

    // First, we fill in the volume-interior nodes. This preserves the original
    // curvature of the mesh.
    const int nElmt = m_element[m_expDim].size();
    int tmpId = 0;
    for (int i = 0; i < nElmt; ++i)
    {
        if (m_verbose)
        {
            LibUtilities::PrintProgressbar(i, nElmt, "MakeOrder: Elements: ");
        }
        ElementSharedPtr el                    = m_element[m_expDim][i];
        SpatialDomains::GeometrySharedPtr geom = el->GetGeom(m_spaceDim);
        geom->FillGeom();
        el->MakeOrder(order, geom, pTypes[el->GetConf().m_e], m_spaceDim, tmpId);
    }

    // Now make copies of each of the edges.
Michael Turner's avatar
Michael Turner committed
164
    for (eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
165
    {
Michael Turner's avatar
Michael Turner committed
166 167
        edgeCopies[(*eit)->m_id] = EdgeSharedPtr(new Edge(*(*eit)));
        (*eit)->m_edgeNodes.clear();
168 169
    }

170 171 172
    // Now copy faces. Make sure that this is a "deep copy", so that the face's
    // edge list corresponds to the copied edges, otherwise we end up in a
    // non-consistent state.
Michael Turner's avatar
Michael Turner committed
173
    for (fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++)
174
    {
175 176 177 178 179 180 181 182
        FaceSharedPtr tmpFace = FaceSharedPtr(new Face(*(*fit)));

        for (int i = 0; i < tmpFace->m_edgeList.size(); ++i)
        {
            tmpFace->m_edgeList[i] = edgeCopies[tmpFace->m_edgeList[i]->m_id];
        }

        faceCopies[(*fit)->m_id] = tmpFace;
Michael Turner's avatar
Michael Turner committed
183
        (*fit)->m_faceNodes.clear();
184 185 186 187
    }

    boost::unordered_set<int> processedEdges, processedFaces, processedVolumes;

188 189
    // note if CAD previously existed on the face or edge, the new points need
    // to be projected onto the CAD entity.
190

191 192
    // Call MakeOrder with our generated geometries on each edge to fill in edge
    // interior nodes.
Michael Turner's avatar
Michael Turner committed
193
    int ct = 0;
194
    for (eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++, ct++)
195
    {
196
        if (m_verbose)
Michael Turner's avatar
Michael Turner committed
197
        {
Michael Turner's avatar
Michael Turner committed
198 199
            LibUtilities::PrintProgressbar(ct, m_edgeSet.size(),
                                           "MakeOrder: Edges: ");
Michael Turner's avatar
Michael Turner committed
200
        }
201 202 203 204 205 206 207
        int edgeId = (*eit)->m_id;

        if (processedEdges.find(edgeId) != processedEdges.end())
        {
            continue;
        }

Michael Turner's avatar
Michael Turner committed
208 209 210 211 212 213
        EdgeSharedPtr cpEdge                   = edgeCopies[edgeId];
        SpatialDomains::GeometrySharedPtr geom = cpEdge->GetGeom(m_spaceDim);
        geom->FillGeom();

        (*eit)->MakeOrder(order, geom, pTypes[LibUtilities::eSegment],
                          m_spaceDim, id);
214 215 216
        processedEdges.insert(edgeId);
    }

217 218
    // Call MakeOrder with our generated geometries on each face to fill in face
    // interior nodes.
Michael Turner's avatar
Michael Turner committed
219
    ct = 0;
220
    for (fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++, ct++)
221
    {
222
        if (m_verbose)
Michael Turner's avatar
Michael Turner committed
223
        {
Michael Turner's avatar
Michael Turner committed
224 225
            LibUtilities::PrintProgressbar(ct, m_faceSet.size(),
                                           "MakeOrder: Faces: ");
Michael Turner's avatar
Michael Turner committed
226
        }
227
        int faceId = (*fit)->m_id;
228

229
        if (processedFaces.find(faceId) != processedFaces.end())
230 231 232 233
        {
            continue;
        }

Michael Turner's avatar
Michael Turner committed
234 235 236 237 238 239 240 241
        FaceSharedPtr cpFace                   = faceCopies[faceId];
        SpatialDomains::GeometrySharedPtr geom = cpFace->GetGeom(m_spaceDim);
        geom->FillGeom();

        LibUtilities::ShapeType type = (*fit)->m_vertexList.size() == 3
                                           ? LibUtilities::eTriangle
                                           : LibUtilities::eQuadrilateral;
        (*fit)->MakeOrder(order, geom, pTypes[type], m_spaceDim, id);
242 243 244
        processedFaces.insert(faceId);
    }

245 246 247 248
    // Copy curvature into boundary conditions
    for (int i = 0; i < m_element[1].size(); ++i)
    {
        ElementSharedPtr el = m_element[1][i];
Michael Turner's avatar
Michael Turner committed
249
        EdgeSharedPtr edge  = el->GetEdgeLink();
250 251 252 253 254 255 256 257 258

        if (!edge)
        {
            continue;
        }

        // Copy face curvature
        el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
                      pTypes[el->GetConf().m_e], m_spaceDim, id, true);
259
        el->SetVolumeNodes(edge->m_edgeNodes);
260 261 262 263 264
    }

    for (int i = 0; i < m_element[2].size(); ++i)
    {
        ElementSharedPtr el = m_element[2][i];
Michael Turner's avatar
Michael Turner committed
265
        FaceSharedPtr face  = el->GetFaceLink();
266 267 268 269 270 271 272 273 274 275 276 277

        if (!face)
        {
            continue;
        }

        // Copy face curvature
        el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
                      pTypes[el->GetConf().m_e], m_spaceDim, id, true);
        el->SetVolumeNodes(face->m_faceNodes);
    }

278 279
    for (int i = 0; i < nElmt; ++i)
    {
280 281
        vector<NodeSharedPtr> tmp = m_element[m_expDim][i]->GetVolumeNodes();
        for (int j = 0; j < tmp.size(); ++j)
Michael Turner's avatar
Michael Turner committed
282
        {
283
            tmp[j]->m_id = id++;
Michael Turner's avatar
Michael Turner committed
284
        }
285
    }
Michael Turner's avatar
Michael Turner committed
286

287 288
    if (m_verbose)
    {
Michael Turner's avatar
Michael Turner committed
289
        cout << endl;
290
    }
291
}
Michael Turner's avatar
Michael Turner committed
292 293
}
}