2DGenerator.cpp 18.4 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.h>
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()
{
75 76
    // Check that cad is 2D
    Array<OneD, NekDouble> bndBox = m_mesh->m_cad->GetBoundingBox();
Julian Marcon's avatar
Julian Marcon committed
77
    ASSERTL0(fabs(bndBox[5] - bndBox[4]) < 1.0e-7, "CAD isn't 2D");
cadfix's avatar
cadfix committed
78

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

Julian Marcon's avatar
Julian Marcon committed
90
    // find the ends of the BL curves
91 92 93 94 95
    if (m_config["blcurves"].beenSet)
    {
        FindBLEnds();
    }

96 97 98 99 100 101 102 103
    // 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");
        }
Michael Turner's avatar
Michael Turner committed
104 105 106
        vector<unsigned int>::iterator f =
            find(m_blCurves.begin(), m_blCurves.end(), i);
        if (f == m_blCurves.end())
107 108 109
        {
            m_curvemeshes[i] =
                MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
110

Julian Marcon's avatar
Julian Marcon committed
111 112 113
            // Fheck if this curve is at an end of the BL
            // If so, define an offset for the second node, corresponding to the
            // BL thickness
114 115 116 117 118 119 120
            if (m_blends.count(i))
            {
                vector<CADVertSharedPtr> vertices =
                    m_mesh->m_cad->GetCurve(i)->GetVertex();
                Array<OneD, NekDouble> loc;
                NekDouble t;

Julian Marcon's avatar
Julian Marcon committed
121
                // offset needed at first node (or both)
Julian Marcon's avatar
Tidy  
Julian Marcon committed
122
                if (m_blends[i] == 0 || m_blends[i] == 2)
123 124 125 126 127 128
                {
                    loc = vertices[0]->GetLoc();
                    t   = m_thickness.Evaluate(m_thickness_ID, loc[0], loc[1],
                                             loc[2], 0.0);
                    m_curvemeshes[i]->SetOffset(0, t);
                }
Julian Marcon's avatar
Julian Marcon committed
129
                // offset needed at second node (or both)
Julian Marcon's avatar
Tidy  
Julian Marcon committed
130
                if (m_blends[i] == 1 || m_blends[i] == 2)
131
                {
Julian Marcon's avatar
Tidy  
Julian Marcon committed
132
                    loc = vertices[1]->GetLoc();
133 134
                    t   = m_thickness.Evaluate(m_thickness_ID, loc[0], loc[1],
                                             loc[2], 0.0);
Julian Marcon's avatar
Tidy  
Julian Marcon committed
135
                    m_curvemeshes[i]->SetOffset(1, t);
136 137
                }
            }
138 139 140
        }
        else
        {
Michael Turner's avatar
Michael Turner committed
141 142
            m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
                i, m_mesh, m_config["blthick"].as<string>());
143
        }
144 145 146
        m_curvemeshes[i]->Mesh();
    }

147 148 149 150 151 152 153 154 155
    ////////
    // consider periodic curves

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

156 157
    ////////////////////////////////////////

158
    if (m_config["blcurves"].beenSet)
Michael Turner's avatar
Michael Turner committed
159
    {
160
        // we need to do the boundary layer generation in a face by face basis
Michael Turner's avatar
Michael Turner committed
161 162 163
        MakeBLPrep();
        for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
        {
Michael Turner's avatar
Michael Turner committed
164
            MakeBL(i);
Michael Turner's avatar
Michael Turner committed
165
        }
166

Julian Marcon's avatar
Julian Marcon committed
167 168
        // If the BL doesn't form closed loops, we need to remove the outside
        // nodes from the curve meshes
169 170 171 172 173 174 175 176 177 178 179 180
        for (map<unsigned, unsigned>::iterator ic = m_blends.begin();
             ic != m_blends.end(); ++ic)
        {
            vector<NodeSharedPtr> nodes =
                m_curvemeshes[ic->first]->GetMeshPoints();

            if (ic->second == 0 || ic->second == 2)
            {
                nodes.erase(nodes.begin());
            }
            if (ic->second == 1 || ic->second == 2)
            {
Julian Marcon's avatar
Julian Marcon committed
181
                nodes.erase(nodes.end() - 1);
182 183
            }

Julian Marcon's avatar
Julian Marcon committed
184 185
            // Rebuild the curvemesh without the first node, the last node or
            // both
186 187 188 189
            m_curvemeshes[ic->first] =
                MemoryManager<CurveMesh>::AllocateSharedPtr(ic->first, m_mesh,
                                                            nodes);
        }
190
    }
Michael Turner's avatar
Michael Turner committed
191
    
192 193 194 195 196
    if (m_mesh->m_verbose)
    {
        cout << endl << "\tFace meshing:" << endl << endl;
    }
    // linear mesh all surfaces
Michael Turner's avatar
Michael Turner committed
197
    for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
198 199 200 201 202 203
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumSurf(),
                                           "Face progress");
        }
204

Julian Marcon's avatar
Julian Marcon committed
205 206
        m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
            i, m_mesh, m_curvemeshes, 99 + i);
Michael Turner's avatar
Michael Turner committed
207
        m_facemeshes[i]->Mesh();
208 209
    }

210 211
    ////////////////////////////////////

212
    EdgeSet::iterator it;
213 214 215 216 217 218 219 220 221 222 223 224 225 226
    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);
227
    }
228

229 230 231 232 233 234 235 236
    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
    ProcessComposites();
    Report();
}

Julian Marcon's avatar
Julian Marcon committed
237 238
void Generator2D::FindBLEnds()
{
Julian Marcon's avatar
Julian Marcon committed
239 240 241 242 243
    // Set of CAD vertices
    // Vertices of each curve are added to the set if not found and removed from
    // the set if found
    // This leaves us with a set of vertices that are at the end of BL open
    // loops
244 245
    set<CADVertSharedPtr> cadverts;

Julian Marcon's avatar
Tidy  
Julian Marcon committed
246
    for (int it = 0; it < m_blCurves.size(); ++it)
Julian Marcon's avatar
Julian Marcon committed
247 248
    {
        vector<CADVertSharedPtr> vertices =
Julian Marcon's avatar
Tidy  
Julian Marcon committed
249
            m_mesh->m_cad->GetCurve(m_blCurves[it])->GetVertex();
Julian Marcon's avatar
Julian Marcon committed
250

Julian Marcon's avatar
Tidy  
Julian Marcon committed
251
        for (int iv = 0; iv < vertices.size(); ++iv)
Julian Marcon's avatar
Julian Marcon committed
252
        {
Julian Marcon's avatar
Tidy  
Julian Marcon committed
253
            set<CADVertSharedPtr>::iterator is = cadverts.find(vertices[iv]);
254 255 256 257 258 259 260

            if (is != cadverts.end())
            {
                cadverts.erase(is);
            }
            else
            {
Julian Marcon's avatar
Tidy  
Julian Marcon committed
261
                cadverts.insert(vertices[iv]);
262 263 264 265
            }
        }
    }

Julian Marcon's avatar
Julian Marcon committed
266 267 268
    // Build m_blends based on the previously constructed set of vertices
    // m_blends is a map of curve number (the curves right outside the BL open
    // loops) to the offset node number: 0, 1 or 2 (for both)
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
    for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); ++i)
    {
        if (find(m_blCurves.begin(), m_blCurves.end(), i) != m_blCurves.end())
        {
            continue;
        }

        vector<CADVertSharedPtr> vertices =
            m_mesh->m_cad->GetCurve(i)->GetVertex();

        for (int j = 0; j < 2; ++j)
        {
            if (!cadverts.count(vertices[j]))
            {
                continue;
            }
Julian Marcon's avatar
Julian Marcon committed
285

286
            if (m_blends.count(i))
Julian Marcon's avatar
Julian Marcon committed
287
            {
288
                m_blends[i] = 2;
Julian Marcon's avatar
Julian Marcon committed
289 290 291
            }
            else
            {
292
                m_blends[i] = j;
Julian Marcon's avatar
Julian Marcon committed
293 294 295 296 297
            }
        }
    }
}

Michael Turner's avatar
Michael Turner committed
298
void Generator2D::MakeBLPrep()
Michael Turner's avatar
Michael Turner committed
299
{
300 301 302 303 304 305 306
    if (m_mesh->m_verbose)
    {
        cout << endl << "\tBoundary layer meshing:" << endl << endl;
    }

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

Julian Marcon's avatar
Tidy  
Julian Marcon committed
307
    for (int it = 0; it < m_blCurves.size(); ++it)
308
    {
Julian Marcon's avatar
Tidy  
Julian Marcon committed
309 310
        vector<EdgeSharedPtr> localedges =
            m_curvemeshes[m_blCurves[it]]->GetMeshEdges();
311
        for (int i = 0; i < localedges.size(); i++)
312
        {
Michael Turner's avatar
Michael Turner committed
313 314
            m_nodesToEdge[localedges[i]->m_n1].push_back(localedges[i]);
            m_nodesToEdge[localedges[i]->m_n2].push_back(localedges[i]);
315 316
        }
    }
Michael Turner's avatar
Michael Turner committed
317
}
318

Michael Turner's avatar
Michael Turner committed
319
void Generator2D::MakeBL(int faceid)
Michael Turner's avatar
Michael Turner committed
320
{
321 322
    map<int, Array<OneD, NekDouble> > edgeNormals;
    int eid = 0;
Julian Marcon's avatar
Tidy  
Julian Marcon committed
323
    for (int it = 0; it < m_blCurves.size(); ++it)
324
    {
325
        CADOrientation::Orientation edgeo =
Julian Marcon's avatar
Tidy  
Julian Marcon committed
326 327 328
            m_mesh->m_cad->GetCurve(m_blCurves[it])->GetOrienationWRT(faceid);
        vector<EdgeSharedPtr> es =
            m_curvemeshes[m_blCurves[it]]->GetMeshEdges();
329 330
        // on each !!!EDGE!!! calculate a normal
        // always to the left unless edgeo is 1
331 332
        // normal must be done in the parametric space (and then projected back)
        // because of face orientation
Michael Turner's avatar
Michael Turner committed
333
        for (int j = 0; j < es.size(); j++)
334 335
        {
            es[j]->m_id = eid++;
Michael Turner's avatar
Michael Turner committed
336
            Array<OneD, NekDouble> p1, p2;
337 338
            p1 = es[j]->m_n1->GetCADSurfInfo(faceid);
            p2 = es[j]->m_n2->GetCADSurfInfo(faceid);
339
            Array<OneD, NekDouble> n(2);
340 341
            n[0] = p1[1] - p2[1];
            n[1] = p2[0] - p1[0];
342
            if (edgeo == CADOrientation::eBackwards)
Michael Turner's avatar
Michael Turner committed
343
            {
Michael Turner's avatar
Michael Turner committed
344 345
                n[0] *= -1.0;
                n[1] *= -1.0;
Michael Turner's avatar
Michael Turner committed
346
            }
Michael Turner's avatar
Michael Turner committed
347
            NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
348 349
            n[0] /= mag;
            n[1] /= mag;
Michael Turner's avatar
Michael Turner committed
350 351 352
            Array<OneD, NekDouble> np(2);
            np[0] = p1[0] + n[0];
            np[1] = p1[1] + n[1];
353 354 355 356 357 358 359
            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;
360 361 362
            edgeNormals[es[j]->m_id] = n;
        }
    }
Julian Marcon's avatar
Julian Marcon committed
363

Michael Turner's avatar
Michael Turner committed
364 365 366
    bool adjust           = m_config["bltadjust"].beenSet;
    NekDouble divider     = m_config["bltadjust"].as<NekDouble>();
    bool adjustEverywhere = m_config["adjustblteverywhere"].beenSet;
367

Michael Turner's avatar
Michael Turner committed
368 369 370 371 372 373
    if (divider < 2.0)
    {
        WARNINGL1(false, "BndLayerAdjustment too low, corrected to 2.0");
        divider = 2.0;
    }

374 375
    map<NodeSharedPtr, NodeSharedPtr> nodeNormals;
    map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
376
    for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
377
    {
378 379 380
        ASSERTL0(it->second.size() == 1 || it->second.size() == 2,
                 "weirdness, most likely bl_surfs are incorrect");

Julian Marcon's avatar
Julian Marcon committed
381 382
        // If node at the end of the BL open loop, the "normal node" isn't
        // constructed by computing a normal but found on the adjacent curve
383 384
        if (it->second.size() == 1)
        {
cadfix's avatar
cadfix committed
385
            vector<CADCurveSharedPtr> curves = it->first->GetCADCurves();
386 387

            vector<EdgeSharedPtr> edges =
cadfix's avatar
cadfix committed
388
                m_curvemeshes[curves[0]->GetId()]->GetMeshEdges();
389 390 391
            vector<EdgeSharedPtr>::iterator ie =
                find(edges.begin(), edges.end(), it->second[0]);
            int rightCurve =
cadfix's avatar
cadfix committed
392
                (ie == edges.end()) ? curves[0]->GetId() : curves[1]->GetId();
393 394 395 396 397 398 399 400 401

            vector<NodeSharedPtr> nodes =
                m_curvemeshes[rightCurve]->GetMeshPoints();
            nodeNormals[it->first] =
                (nodes[0] == it->first) ? nodes[1] : nodes[nodes.size() - 2];

            continue;
        }

402
        Array<OneD, NekDouble> n(3, 0.0);
403 404
        Array<OneD, NekDouble> n1 = edgeNormals[it->second[0]->m_id];
        Array<OneD, NekDouble> n2 = edgeNormals[it->second[1]->m_id];
Michael Turner's avatar
Michael Turner committed
405 406 407
        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]);
408 409
        n[0] /= mag;
        n[1] /= mag;
Michael Turner's avatar
Michael Turner committed
410 411
        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
412 413 414 415 416 417 418 419 420 421
        // 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);
            }
        }
Michael Turner's avatar
Michael Turner committed
422

423 424
        n[0]             = n[0] * t + it->first->m_x;
        n[1]             = n[1] * t + it->first->m_y;
425
        NodeSharedPtr nn = std::shared_ptr<Node>(
Dave Moxey's avatar
Dave Moxey committed
426
            new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
427
        CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
mt4313's avatar
mt4313 committed
428
        Array<OneD, NekDouble> uv = s->locuv(n);
429
        nn->SetCADSurf(s, uv);
430 431
        nodeNormals[it->first] = nn;
    }
432

Julian Marcon's avatar
Tidy  
Julian Marcon committed
433
    for (int it = 0; it < m_blCurves.size(); ++it)
434
    {
435
        CADOrientation::Orientation edgeo =
Julian Marcon's avatar
Tidy  
Julian Marcon committed
436 437 438
            m_mesh->m_cad->GetCurve(m_blCurves[it])->GetOrienationWRT(faceid);
        vector<NodeSharedPtr> ns =
            m_curvemeshes[m_blCurves[it]]->GetMeshPoints();
439
        vector<NodeSharedPtr> newNs;
Michael Turner's avatar
Michael Turner committed
440
        for (int i = 0; i < ns.size(); i++)
441 442 443
        {
            newNs.push_back(nodeNormals[ns[i]]);
        }
Julian Marcon's avatar
Tidy  
Julian Marcon committed
444 445 446
        m_curvemeshes[m_blCurves[it]] =
            MemoryManager<CurveMesh>::AllocateSharedPtr(m_blCurves[it], m_mesh,
                                                        newNs);
447
        if (edgeo == CADOrientation::eBackwards)
448 449 450 451 452 453 454
        {
            reverse(ns.begin(), ns.end());
        }
        for (int i = 0; i < ns.size() - 1; ++i)
        {
            vector<NodeSharedPtr> qns;
            qns.push_back(ns[i]);
Michael Turner's avatar
Michael Turner committed
455 456
            qns.push_back(ns[i + 1]);
            qns.push_back(nodeNormals[ns[i + 1]]);
457 458 459
            qns.push_back(nodeNormals[ns[i]]);
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
            vector<int> tags;
460
            tags.push_back(101);
461 462 463
            ElementSharedPtr E = GetElementFactory().CreateInstance(
                LibUtilities::eQuadrilateral, conf, qns, tags);
            E->m_parentCAD = m_mesh->m_cad->GetSurf(faceid);
464 465
            for (int j = 0; j < E->GetEdgeCount(); ++j)
            {
Michael Turner's avatar
Michael Turner committed
466
                pair<EdgeSet::iterator, bool> testIns;
467 468 469 470 471 472 473 474 475
                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);
                }
            }
476
            m_mesh->m_element[2].push_back(E);
Julian Marcon's avatar
Julian Marcon committed
477
        }
478
    }
Michael Turner's avatar
Michael Turner committed
479 480
}

481 482 483 484 485 486 487 488 489 490 491
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(":"));

Julian Marcon's avatar
Tidy  
Julian Marcon committed
492
    for (int il = 0; il < lines.size(); ++il)
493 494
    {
        vector<string> tmp;
Julian Marcon's avatar
Tidy  
Julian Marcon committed
495
        boost::split(tmp, lines[il], boost::is_any_of(","));
496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518

        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
519
        m_curvemeshes[ip->second]->PeriodicOverwrite(m_curvemeshes[ip->first]);
520 521 522 523 524 525 526 527 528 529 530 531 532 533 534
    }

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

535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
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;
    }
}
550 551
}
}