2DGenerator.cpp 18.5 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()
{
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");
78
    
Michael Turner's avatar
Michael Turner committed
79
    if (m_mesh->m_verbose)
80
    {
Michael Turner's avatar
Michael Turner committed
81
82
83
84
85
86
87
88
        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);
89

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

Michael Turner's avatar
Michael Turner committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    // 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);
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
                }
            }
Michael Turner's avatar
Michael Turner committed
138
139
140
141
142
143
144
        }
        else
        {
            m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
                i, m_mesh, m_config["blthick"].as<string>());
        }
        m_curvemeshes[i]->Mesh();
145
146
    }

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

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

Michael Turner's avatar
Michael Turner committed
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
191
    }

192
    if (m_mesh->m_verbose)
193
194
195
196
    {
        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
    }
Michael Turner's avatar
Michael Turner committed
209
210
211

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

212
    EdgeSet::iterator it;
Michael Turner's avatar
Michael Turner committed
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
    }
Michael Turner's avatar
Michael Turner committed
228

229
230
    // m_mesh->m_expDim = 1;
    // m_mesh->m_element[2].clear();
Michael Turner's avatar
Michael Turner committed
231

232
233
234
235
236
237
238
239
    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
    ProcessComposites();
    Report();
}

Julian Marcon's avatar
Julian Marcon committed
240
241
void Generator2D::FindBLEnds()
{
Julian Marcon's avatar
Julian Marcon committed
242
243
244
245
246
    // 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
247
248
    set<CADVertSharedPtr> cadverts;

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

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

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

Julian Marcon's avatar
Julian Marcon committed
269
270
271
    // 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)
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    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
288

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

Michael Turner's avatar
Michael Turner committed
301
void Generator2D::MakeBLPrep()
Michael Turner's avatar
Michael Turner committed
302
{
303
304
305
306
307
308
309
    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
310
    for (int it = 0; it < m_blCurves.size(); ++it)
311
    {
Julian Marcon's avatar
Tidy    
Julian Marcon committed
312
313
        vector<EdgeSharedPtr> localedges =
            m_curvemeshes[m_blCurves[it]]->GetMeshEdges();
314
        for (int i = 0; i < localedges.size(); i++)
315
        {
Michael Turner's avatar
Michael Turner committed
316
317
            m_nodesToEdge[localedges[i]->m_n1].push_back(localedges[i]);
            m_nodesToEdge[localedges[i]->m_n2].push_back(localedges[i]);
318
319
        }
    }
Michael Turner's avatar
Michael Turner committed
320
}
321

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

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

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

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

Julian Marcon's avatar
Julian Marcon committed
384
385
        // 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
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
        if (it->second.size() == 1)
        {
            vector<pair<int, CADCurveSharedPtr> > curves =
                it->first->GetCADCurves();

            vector<EdgeSharedPtr> edges =
                m_curvemeshes[curves[0].first]->GetMeshEdges();
            vector<EdgeSharedPtr>::iterator ie =
                find(edges.begin(), edges.end(), it->second[0]);
            int rightCurve =
                (ie == edges.end()) ? curves[0].first : curves[1].first;

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

            continue;
        }

406
        Array<OneD, NekDouble> n(3, 0.0);
407
408
        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
409
410
411
        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]);
412
413
        n[0] /= mag;
        n[1] /= mag;
Michael Turner's avatar
Michael Turner committed
414
415
        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
416
417
418
419
420
421
422
423
424
425
        // 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
426

427
428
        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
429
430
        NodeSharedPtr nn = boost::shared_ptr<Node>(
            new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
431
432
        CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
        Array<OneD, NekDouble> uv = s->locuv(n);
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
433
        nn->SetCADSurf(faceid, s, uv);
434
435
        nodeNormals[it->first] = nn;
    }
436

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

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

        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
523
        m_curvemeshes[ip->second]->PeriodicOverwrite(m_curvemeshes[ip->first]);
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
    }

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

539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
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;
    }
}
554
555
}
}