Mesh.cpp 9.18 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
{
Michael Turner's avatar
Michael Turner committed
101 102 103 104 105 106 107

    // going to make a copy of the curavture information
    // cheaper that geom 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

108 109 110 111 112
    int id = m_vertexSet.size();

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

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

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

Michael Turner's avatar
Michael Turner committed
144 145
    // Begin by coping mesh objects for edges, faces
    // so that we don't affect any neighbouring elements in the mesh as
146
    // we process each element.
Michael Turner's avatar
Michael Turner committed
147 148
    // at the same time we delete the curvature from the original edge
    for (eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
149
    {
Michael Turner's avatar
Michael Turner committed
150 151
        edgeCopies[(*eit)->m_id] = EdgeSharedPtr(new Edge(*(*eit)));
        (*eit)->m_edgeNodes.clear();
152 153
    }

Michael Turner's avatar
Michael Turner committed
154
    for (fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++)
155
    {
Michael Turner's avatar
Michael Turner committed
156 157
        faceCopies[(*fit)->m_id] = FaceSharedPtr(new Face(*(*fit)));
        (*fit)->m_faceNodes.clear();
158 159 160 161
    }

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

162 163
    // note if CAD previously existed on the face or edge, the new points need
    // to be projected onto the CAD entity.
164

165 166
    // Call MakeOrder with our generated geometries on each edge to fill in edge
    // interior nodes.
Michael Turner's avatar
Michael Turner committed
167
    int ct = 0;
168
    for (eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++, ct++)
169
    {
170
        if (m_verbose)
Michael Turner's avatar
Michael Turner committed
171
        {
Michael Turner's avatar
Michael Turner committed
172 173
            LibUtilities::PrintProgressbar(ct, m_edgeSet.size(),
                                           "MakeOrder: Edges: ");
Michael Turner's avatar
Michael Turner committed
174
        }
175 176 177 178 179 180 181
        int edgeId = (*eit)->m_id;

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

Michael Turner's avatar
Michael Turner committed
182 183 184 185 186 187
        EdgeSharedPtr cpEdge                   = edgeCopies[edgeId];
        SpatialDomains::GeometrySharedPtr geom = cpEdge->GetGeom(m_spaceDim);
        geom->FillGeom();

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

191 192
    // Call MakeOrder with our generated geometries on each face to fill in face
    // interior nodes.
Michael Turner's avatar
Michael Turner committed
193
    ct = 0;
194
    for (fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++, 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_faceSet.size(),
                                           "MakeOrder: Faces: ");
Michael Turner's avatar
Michael Turner committed
200
        }
201
        int faceId = (*fit)->m_id;
202

203
        if (processedFaces.find(faceId) != processedFaces.end())
204 205 206 207
        {
            continue;
        }

Michael Turner's avatar
Michael Turner committed
208 209 210 211 212 213 214 215
        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);
216 217 218
        processedFaces.insert(faceId);
    }

219 220 221 222
    // 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
223
        EdgeSharedPtr edge  = el->GetEdgeLink();
224 225 226 227 228 229 230 231 232

        if (!edge)
        {
            continue;
        }

        // Copy face curvature
        el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
                      pTypes[el->GetConf().m_e], m_spaceDim, id, true);
233
        el->SetVolumeNodes(edge->m_edgeNodes);
234 235 236 237 238
    }

    for (int i = 0; i < m_element[2].size(); ++i)
    {
        ElementSharedPtr el = m_element[2][i];
Michael Turner's avatar
Michael Turner committed
239
        FaceSharedPtr face  = el->GetFaceLink();
240 241 242 243 244 245 246 247 248 249 250 251 252

        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);
    }

    // Finally, fill in volumes.
253 254 255
    const int nElmt = m_element[m_expDim].size();
    for (int i = 0; i < nElmt; ++i)
    {
256
        if (m_verbose)
Michael Turner's avatar
Michael Turner committed
257
        {
258
            LibUtilities::PrintProgressbar(i, nElmt, "MakeOrder: Elements: ");
Michael Turner's avatar
Michael Turner committed
259
        }
Michael Turner's avatar
Michael Turner committed
260 261 262 263
        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, id);
264
    }
Michael Turner's avatar
Michael Turner committed
265

266 267
    if (m_verbose)
    {
Michael Turner's avatar
Michael Turner committed
268
        cout << endl;
269
    }
270
}
Michael Turner's avatar
Michael Turner committed
271 272
}
}