2DGenerator.cpp 13.7 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
//  File: Generator2D.cpp
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.
//
32
//  Description: 2D generator object methods.
33 34 35
//
////////////////////////////////////////////////////////////////////////////////
#include <algorithm>
Marcon's avatar
Marcon committed
36
#include <math.h>
37 38 39

#include <NekMeshUtils/2DGenerator/2DGenerator.h>

40
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
41
#include <LibUtilities/BasicUtils/Progressbar.hpp>
42

43 44
#include <boost/algorithm/string.hpp>

45 46 47 48 49 50 51
using namespace std;
namespace Nektar
{
namespace NekMeshUtils
{

ModuleKey Generator2D::className = GetModuleFactory().RegisterCreatorFunction(
52
    ModuleKey(eProcessModule, "2dgenerator"), Generator2D::create,
53 54 55 56
    "Generates a 2D mesh");

Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
{
57
    m_config["blcurves"] =
Michael Turner's avatar
Michael Turner committed
58
        ConfigOption(false, "", "Generate parallelograms on these curves");
59
    m_config["blthick"] =
Michael Turner's avatar
Michael Turner committed
60
        ConfigOption(false, "0.0", "Parallelogram layer thickness");
61
    m_config["periodic"] =
Michael Turner's avatar
Michael Turner committed
62
        ConfigOption(false, "", "Set of pairs of periodic curves");
63 64 65 66
    m_config["bltadjust"] =
        ConfigOption(false, "2.0", "Boundary layer thickness adjustment");
    m_config["adjustblteverywhere"] =
        ConfigOption(true, "0", "Adjust thickness everywhere");
67 68 69 70 71 72 73 74
}

Generator2D::~Generator2D()
{
}

void Generator2D::Process()
{
Michael Turner's avatar
Michael Turner committed
75
    if (m_mesh->m_verbose)
76
    {
Michael Turner's avatar
Michael Turner committed
77 78 79 80 81 82 83 84
        cout << endl << "2D meshing" << endl;
        cout << endl << "\tCurve meshing:" << endl << endl;
    }
    m_mesh->m_numNodes = m_mesh->m_cad->GetNumVerts();
    m_thickness_ID =
        m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
    ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
                                  m_blCurves);
85

Michael Turner's avatar
Michael Turner committed
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
    // linear mesh all curves
    for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumCurve(),
                                           "Curve progress");
        }
        vector<unsigned int>::iterator f =
            find(m_blCurves.begin(), m_blCurves.end(), i);
        if (f == m_blCurves.end())
        {
            m_curvemeshes[i] =
                MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
        }
        else
        {
            m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
                i, m_mesh, m_config["blthick"].as<string>());
        }
        m_curvemeshes[i]->Mesh();
107 108
    }

109 110 111 112 113 114 115 116 117
    ////////
    // consider periodic curves

    if (m_config["periodic"].beenSet)
    {
        PeriodicPrep();
        MakePeriodic();
    }

Michael Turner's avatar
Michael Turner committed
118 119
    ////////////////////////////////////////

120
    if (m_config["blcurves"].beenSet)
Michael Turner's avatar
Michael Turner committed
121
    {
122
        // we need to do the boundary layer generation in a face by face basis
Michael Turner's avatar
Michael Turner committed
123 124 125
        MakeBLPrep();
        for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
        {
Michael Turner's avatar
Michael Turner committed
126
            MakeBL(i);
Michael Turner's avatar
Michael Turner committed
127
        }
128 129
    }

130
    if (m_mesh->m_verbose)
131 132 133 134
    {
        cout << endl << "\tFace meshing:" << endl << endl;
    }
    // linear mesh all surfaces
Michael Turner's avatar
Michael Turner committed
135
    for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
136 137 138 139 140 141
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumSurf(),
                                           "Face progress");
        }
142

Julian Marcon's avatar
Julian Marcon committed
143 144
        m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
            i, m_mesh, m_curvemeshes, 99 + i);
Michael Turner's avatar
Michael Turner committed
145
        m_facemeshes[i]->Mesh();
146
    }
Michael Turner's avatar
Michael Turner committed
147 148 149

    ////////////////////////////////////

150
    EdgeSet::iterator it;
Michael Turner's avatar
Michael Turner committed
151 152 153 154 155 156 157 158 159 160 161 162 163 164
    for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
    {
        vector<NodeSharedPtr> ns;
        ns.push_back((*it)->m_n1);
        ns.push_back((*it)->m_n2);
        // for each iterator create a LibUtilities::eSegement
        // push segment into m_mesh->m_element[1]
        // tag for the elements shoudl be the CAD number of the curves
        ElmtConfig conf(LibUtilities::eSegment, 1, false, false);
        vector<int> tags;
        tags.push_back((*it)->m_parentCAD->GetId());
        ElementSharedPtr E2 = GetElementFactory().CreateInstance(
            LibUtilities::eSegment, conf, ns, tags);
        m_mesh->m_element[1].push_back(E2);
165
    }
Michael Turner's avatar
Michael Turner committed
166

167 168
    // m_mesh->m_expDim = 1;
    // m_mesh->m_element[2].clear();
Michael Turner's avatar
Michael Turner committed
169

170 171 172 173 174 175 176 177
    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
    ProcessComposites();
    Report();
}

Michael Turner's avatar
Michael Turner committed
178
void Generator2D::MakeBLPrep()
Michael Turner's avatar
Michael Turner committed
179
{
180 181 182 183 184 185 186
    if (m_mesh->m_verbose)
    {
        cout << endl << "\tBoundary layer meshing:" << endl << endl;
    }

    // identify the nodes which will become the boundary layer.

187 188
    for (vector<unsigned>::iterator it = m_blCurves.begin();
         it != m_blCurves.end(); ++it)
189
    {
190 191
        vector<EdgeSharedPtr> localedges = m_curvemeshes[*it]->GetMeshEdges();
        for (int i = 0; i < localedges.size(); i++)
192
        {
Michael Turner's avatar
Michael Turner committed
193 194
            m_nodesToEdge[localedges[i]->m_n1].push_back(localedges[i]);
            m_nodesToEdge[localedges[i]->m_n2].push_back(localedges[i]);
195 196
        }
    }
Michael Turner's avatar
Michael Turner committed
197
}
198

Michael Turner's avatar
Michael Turner committed
199
void Generator2D::MakeBL(int faceid)
Michael Turner's avatar
Michael Turner committed
200
{
201 202 203
    map<int, Array<OneD, NekDouble> > edgeNormals;
    int eid = 0;
    for (vector<unsigned>::iterator it = m_blCurves.begin();
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
204
         it != m_blCurves.end(); ++it)
205
    {
206
        CADOrientation::Orientation edgeo =
207
            m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
208 209 210
        vector<EdgeSharedPtr> es = m_curvemeshes[*it]->GetMeshEdges();
        // on each !!!EDGE!!! calculate a normal
        // always to the left unless edgeo is 1
211 212
        // normal must be done in the parametric space (and then projected back)
        // because of face orientation
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
213
        for (int j = 0; j < es.size(); j++)
214 215
        {
            es[j]->m_id = eid++;
Michael Turner's avatar
Michael Turner committed
216
            Array<OneD, NekDouble> p1, p2;
217 218
            p1 = es[j]->m_n1->GetCADSurfInfo(faceid);
            p2 = es[j]->m_n2->GetCADSurfInfo(faceid);
219
            Array<OneD, NekDouble> n(2);
220 221
            n[0] = p1[1] - p2[1];
            n[1] = p2[0] - p1[0];
Michael Turner's avatar
Michael Turner committed
222 223 224 225 226
            if (edgeo == CADOrientation::eBackwards)
            {
                n[0] *= -1.0;
                n[1] *= -1.0;
            }
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
227
            NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
228 229
            n[0] /= mag;
            n[1] /= mag;
Michael Turner's avatar
Michael Turner committed
230 231 232
            Array<OneD, NekDouble> np(2);
            np[0] = p1[0] + n[0];
            np[1] = p1[1] + n[1];
233 234 235 236 237 238 239
            Array<OneD, NekDouble> loc  = es[j]->m_n1->GetLoc();
            Array<OneD, NekDouble> locp = m_mesh->m_cad->GetSurf(faceid)->P(np);
            n[0] = locp[0] - loc[0];
            n[1] = locp[1] - loc[1];
            mag  = sqrt(n[0] * n[0] + n[1] * n[1]);
            n[0] /= mag;
            n[1] /= mag;
240 241 242
            edgeNormals[es[j]->m_id] = n;
        }
    }
Julian Marcon's avatar
Julian Marcon committed
243

Michael Turner's avatar
Michael Turner committed
244 245 246
    bool adjust           = m_config["bltadjust"].beenSet;
    NekDouble divider     = m_config["bltadjust"].as<NekDouble>();
    bool adjustEverywhere = m_config["adjustblteverywhere"].beenSet;
247

Michael Turner's avatar
Michael Turner committed
248 249 250 251 252 253
    if (divider < 2.0)
    {
        WARNINGL1(false, "BndLayerAdjustment too low, corrected to 2.0");
        divider = 2.0;
    }

254 255
    map<NodeSharedPtr, NodeSharedPtr> nodeNormals;
    map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
256
    for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
257
    {
258
        Array<OneD, NekDouble> n(3, 0.0);
Michael Turner's avatar
Michael Turner committed
259 260
        ASSERTL0(it->second.size() == 2,
                 "wierdness, most likely bl_surfs are incorrect");
261 262
        Array<OneD, NekDouble> n1 = edgeNormals[it->second[0]->m_id];
        Array<OneD, NekDouble> n2 = edgeNormals[it->second[1]->m_id];
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
263 264 265
        n[0]          = (n1[0] + n2[0]) / 2.0;
        n[1]          = (n1[1] + n2[1]) / 2.0;
        NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
266 267
        n[0] /= mag;
        n[1] /= mag;
Michael Turner's avatar
Michael Turner committed
268 269
        NekDouble t = m_thickness.Evaluate(m_thickness_ID, it->first->m_x,
                                           it->first->m_y, 0.0, 0.0);
Michael Turner's avatar
Michael Turner committed
270 271 272 273 274 275 276 277 278 279
        // Adjust thickness according to angle between normals
        if (adjust)
        {
            if (adjustEverywhere || it->first->GetNumCadCurve() > 1)
            {
                NekDouble angle = acos(n1[0] * n2[0] + n1[1] * n2[1]);
                angle           = (angle > M_PI) ? 2 * M_PI - angle : angle;
                t /= cos(angle / divider);
            }
        }
Julian Marcon's avatar
Julian Marcon committed
280

281 282
        n[0]             = n[0] * t + it->first->m_x;
        n[1]             = n[1] * t + it->first->m_y;
Dave Moxey's avatar
Dave Moxey committed
283 284
        NodeSharedPtr nn = boost::shared_ptr<Node>(
            new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
285 286
        CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
        Array<OneD, NekDouble> uv = s->locuv(n);
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
287
        nn->SetCADSurf(faceid, s, uv);
288 289
        nodeNormals[it->first] = nn;
    }
290

291
    for (vector<unsigned>::iterator it = m_blCurves.begin();
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
292
         it != m_blCurves.end(); ++it)
293
    {
294
        CADOrientation::Orientation edgeo =
295
            m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
296 297
        vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
        vector<NodeSharedPtr> newNs;
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
298
        for (int i = 0; i < ns.size(); i++)
299 300 301
        {
            newNs.push_back(nodeNormals[ns[i]]);
        }
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
302 303
        m_curvemeshes[*it] =
            MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);
304
        if (edgeo == CADOrientation::eBackwards)
305 306 307 308 309 310 311
        {
            reverse(ns.begin(), ns.end());
        }
        for (int i = 0; i < ns.size() - 1; ++i)
        {
            vector<NodeSharedPtr> qns;
            qns.push_back(ns[i]);
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
312 313
            qns.push_back(ns[i + 1]);
            qns.push_back(nodeNormals[ns[i + 1]]);
314 315 316
            qns.push_back(nodeNormals[ns[i]]);
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
            vector<int> tags;
317
            tags.push_back(101);
318 319 320
            ElementSharedPtr E = GetElementFactory().CreateInstance(
                LibUtilities::eQuadrilateral, conf, qns, tags);
            E->m_parentCAD = m_mesh->m_cad->GetSurf(faceid);
321 322
            for (int j = 0; j < E->GetEdgeCount(); ++j)
            {
Julian Marcon's avatar
Tidy.  
Julian Marcon committed
323
                pair<EdgeSet::iterator, bool> testIns;
324 325 326 327 328 329 330 331 332
                EdgeSharedPtr ed = E->GetEdge(j);
                // look for edge in m_mesh edgeset from curves
                EdgeSet::iterator s = m_mesh->m_edgeSet.find(ed);
                if (!(s == m_mesh->m_edgeSet.end()))
                {
                    ed = *s;
                    E->SetEdge(j, *s);
                }
            }
333
            m_mesh->m_element[2].push_back(E);
Julian Marcon's avatar
Julian Marcon committed
334
        }
335
    }
Michael Turner's avatar
Michael Turner committed
336
}
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

void Generator2D::PeriodicPrep()
{
    m_periodicPairs.clear();
    set<unsigned> periodic;

    // Build periodic curve pairs
    string s = m_config["periodic"].as<string>();
    vector<string> lines;

    boost::split(lines, s, boost::is_any_of(":"));

    for (vector<string>::iterator il = lines.begin(); il != lines.end(); ++il)
    {
        vector<string> tmp;
        boost::split(tmp, *il, boost::is_any_of(","));

        ASSERTL0(tmp.size() == 2, "periodic pairs ill-defined");

        vector<unsigned> data(2);
        data[0] = boost::lexical_cast<unsigned>(tmp[0]);
        data[1] = boost::lexical_cast<unsigned>(tmp[1]);

        ASSERTL0(!periodic.count(data[0]), "curve already periodic");
        ASSERTL0(!periodic.count(data[1]), "curve already periodic");

        m_periodicPairs[data[0]] = data[1];
        periodic.insert(data[0]);
        periodic.insert(data[1]);
    }
}

void Generator2D::MakePeriodic()
{
    // Override slave curves

    for (map<unsigned, unsigned>::iterator ip = m_periodicPairs.begin();
         ip != m_periodicPairs.end(); ++ip)
    {
Michael Turner's avatar
Michael Turner committed
376
        m_curvemeshes[ip->second]->PeriodicOverwrite(m_curvemeshes[ip->first]);
377 378 379 380 381 382 383 384 385 386 387 388 389 390 391
    }

    if (m_mesh->m_verbose)
    {
        cout << "\t\tPeriodic boundary conditions" << endl;
        for (map<unsigned, unsigned>::iterator it = m_periodicPairs.begin();
             it != m_periodicPairs.end(); ++it)
        {
            cout << "\t\t\tCurves " << it->first << " => " << it->second
                 << endl;
        }
        cout << endl;
    }
}

392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
void Generator2D::Report()
{
    if (m_mesh->m_verbose)
    {
        int ns = m_mesh->m_vertexSet.size();
        int es = m_mesh->m_edgeSet.size();
        int ts = m_mesh->m_element[2].size();
        int ep = ns - es + ts;
        cout << endl << "\tSurface mesh statistics" << endl;
        cout << "\t\tNodes: " << ns << endl;
        cout << "\t\tEdges: " << es << endl;
        cout << "\t\tTriangles " << ts << endl;
        cout << "\t\tEuler-Poincaré characteristic: " << ep << endl;
    }
}
407 408
}
}