2DGenerator.cpp 16.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
36
37
38
//
////////////////////////////////////////////////////////////////////////////////
#include <algorithm>

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

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

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

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

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

Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
{
56
57
58
59
    m_config["blcurves"] =
        ConfigOption(false, "0", "Generate parallelograms on these curves");
    m_config["blthick"] =
        ConfigOption(false, "0", "Parallelogram layer thickness");
60
61
    m_config["periodic"] =
        ConfigOption(false, "0", "Set of pairs of periodic curves");
62
63
64
65
66
67
68
69
}

Generator2D::~Generator2D()
{
}

void Generator2D::Process()
{
70
71
72
73
74
75
76
77
    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();

78
    set<unsigned> periodic;
79

80
81
82
    if (m_config["periodic"].beenSet)
    {
        m_periodicPairs.clear();
83

Julian Marcon's avatar
Julian Marcon committed
84
        // Build periodic curve pairs
85

86
87
88
89
90
        string s = m_config["periodic"].as<string>();
        vector<string> lines;

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

Julian Marcon's avatar
Julian Marcon committed
91
92
        for (vector<string>::iterator il = lines.begin(); il != lines.end();
             ++il)
93
94
        {
            vector<unsigned> data;
Julian Marcon's avatar
Julian Marcon committed
95
            ParseUtils::GenerateOrderedVector(il->c_str(), data);
96
97

            ASSERTL0(data.size() == 2, "periodic pairs ill-defined");
98
99
            ASSERTL0(!periodic.count(data[0]), "curve already periodic");
            ASSERTL0(!periodic.count(data[1]), "curve already periodic");
100
101

            m_periodicPairs[data[0]] = data[1];
102
103
            periodic.insert(data[0]);
            periodic.insert(data[1]);
104
105
106
107
108
109
110
        }

        // Check compatibility

        for (map<unsigned, unsigned>::iterator it = m_periodicPairs.begin();
             it != m_periodicPairs.end(); ++it)
        {
Julian Marcon's avatar
Julian Marcon committed
111
112
            NekDouble L1 = m_mesh->m_cad->GetCurve(it->first)->GetTotLength();
            NekDouble L2 = m_mesh->m_cad->GetCurve(it->second)->GetTotLength();
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
113

Julian Marcon's avatar
Julian Marcon committed
114
            ASSERTL0(abs((L1 - L2) / L1) < 1.0e-3,
115
116
                     "periodic curves of different length");
        }
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    }

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

        m_curvemeshes[i] =
            MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);

        m_curvemeshes[i]->Mesh();
    }

    if (m_config["periodic"].beenSet)
    {
        // Override slave curves

        for (map<unsigned, unsigned>::iterator ip = m_periodicPairs.begin();
             ip != m_periodicPairs.end(); ++ip)
        {
            Array<OneD, NekDouble> A1 =
                m_curvemeshes[ip->first]->GetFirstPoint()->GetLoc();
            Array<OneD, NekDouble> A2 =
                m_curvemeshes[ip->first]->GetLastPoint()->GetLoc();
            Array<OneD, NekDouble> B1 =
                m_curvemeshes[ip->second]->GetFirstPoint()->GetLoc();
            Array<OneD, NekDouble> B2 =
                m_curvemeshes[ip->second]->GetLastPoint()->GetLoc();

            Array<OneD, NekDouble> T1(2);
            Array<OneD, NekDouble> T2(2);
Julian Marcon's avatar
Julian Marcon committed
152
            Array<OneD, NekDouble> dT(2);
153

Julian Marcon's avatar
Julian Marcon committed
154
155
            // Compute translation vector

156
157
            T1[0] = B1[0] - A1[0];
            T1[1] = B1[1] - A1[1];
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
158

159
160
161
            T2[0] = B2[0] - A2[0];
            T2[1] = B2[1] - A2[1];

Julian Marcon's avatar
Julian Marcon committed
162
163
164
165
166
167
            dT[0] = T2[0] - T1[0];
            dT[1] = T2[1] - T1[1];

            NekDouble dTmag = (dT[0] * dT[0] + dT[1] * dT[1]) /
                              (T1[0] * T1[0] + T1[1] * T1[1]);

Julian Marcon's avatar
Julian Marcon committed
168
169
            // Check if slave vector is reverse oriented

170
171
            bool reverse = false;

Julian Marcon's avatar
Julian Marcon committed
172
            if (dTmag > 1.0e-3)
173
174
175
176
177
            {
                reverse = true;

                T1[0] = B1[0] - A2[0];
                T1[1] = B1[1] - A2[1];
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
178

179
180
181
                T2[0] = B2[0] - A1[0];
                T2[1] = B2[1] - A1[1];

Julian Marcon's avatar
Julian Marcon committed
182
183
184
185
186
187
188
                dT[0] = T2[0] - T1[0];
                dT[1] = T2[1] - T1[1];

                dTmag = (dT[0] * dT[0] + dT[1] * dT[1]) /
                        (T1[0] * T1[0] + T1[1] * T1[1]);

                ASSERTL0(dTmag < 1.0e-3, "curve cannot be translated");
189
190
            }

Julian Marcon's avatar
Julian Marcon committed
191
192
            // Build vector of translated nodes

193
194
195
196
            vector<NodeSharedPtr> nodes =
                m_curvemeshes[ip->first]->GetMeshPoints();
            vector<NodeSharedPtr> nnodes;

Julian Marcon's avatar
Julian Marcon committed
197
198
            vector<CADSurfSharedPtr> surfs =
                m_curvemeshes[ip->second]->GetCADCurve()->GetAdjSurf();
199
200
201
202
203
204
205
206

            nnodes.push_back(m_curvemeshes[ip->second]->GetFirstPoint());

            for (vector<NodeSharedPtr>::iterator in = nodes.begin() + 1;
                 in != nodes.end() - 1; ++in)
            {
                Array<OneD, NekDouble> loc = (*in)->GetLoc();
                NodeSharedPtr nn = boost::shared_ptr<Node>(new Node(
Julian Marcon's avatar
Julian Marcon committed
207
                    m_mesh->m_numNodes++, loc[0] + T1[0], loc[1] + T1[1], 0.0));
208

Julian Marcon's avatar
Julian Marcon committed
209
210
                for (vector<CADSurfSharedPtr>::iterator is = surfs.begin();
                     is != surfs.end(); ++is)
211
                {
Julian Marcon's avatar
Julian Marcon committed
212
213
                    nn->SetCADSurf((*is)->GetId(), *is,
                                   (*is)->locuv(nn->GetLoc()));
214
215
                }

Julian Marcon's avatar
Julian Marcon committed
216
217
218
219
220
                nn->SetCADCurve(ip->second,
                                m_curvemeshes[ip->second]->GetCADCurve(),
                                m_curvemeshes[ip->second]->GetCADCurve()->loct(
                                    nn->GetLoc()));

221
222
223
224
225
                nnodes.push_back(nn);
            }

            nnodes.push_back(m_curvemeshes[ip->second]->GetLastPoint());

Julian Marcon's avatar
Julian Marcon committed
226
227
            // Reverse internal nodes of the vector if necessary

228
229
230
231
            if (reverse)
            {
                vector<NodeSharedPtr> tmp;

Julian Marcon's avatar
Julian Marcon committed
232
                tmp.push_back(nnodes.front());
233
234
235
236
237
238
239
240
241

                for (vector<NodeSharedPtr>::reverse_iterator rin =
                         nnodes.rbegin() + 1;
                     rin != nnodes.rend() - 1; ++rin)
                {
                    tmp.push_back(*rin);
                }

                tmp.push_back(nnodes.back());
Julian Marcon's avatar
Julian Marcon committed
242
243
244
245

                nnodes.swap(tmp);
            }

Julian Marcon's avatar
Julian Marcon committed
246
247
            // Clean m_edgeSet and build new CurveMesh

Julian Marcon's avatar
Julian Marcon committed
248
249
250
251
252
253
            vector<EdgeSharedPtr> edges =
                m_curvemeshes[ip->second]->GetMeshEdges();
            for (vector<EdgeSharedPtr>::iterator ie = edges.begin();
                 ie != edges.end(); ++ie)
            {
                m_mesh->m_edgeSet.erase(*ie);
254
255
256
257
258
259
            }

            m_curvemeshes[ip->second] =
                MemoryManager<CurveMesh>::AllocateSharedPtr(ip->second, m_mesh,
                                                            nnodes);
        }
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
260

261
262
263
264
265
266
        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)
            {
267
268
                cout << "\t\t\tCurves " << it->first << " => " << it->second
                     << endl;
269
270
271
272
273
            }
            cout << endl;
        }
    }

274
275
276
277
278
279
280
281
282
283
284
285
286
287
    EdgeSet::iterator it;
    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;
Michael Turner's avatar
Michael Turner committed
288
        tags.push_back((*it)->m_parentCAD->GetId());
289
290
291
292

        ElementSharedPtr E2 = GetElementFactory().CreateInstance(
            LibUtilities::eSegment, conf, ns, tags);

Michael Turner's avatar
Michael Turner committed
293
        m_mesh->m_element[1].push_back(E2);
294
295
    }

296
    if (m_config["blcurves"].beenSet)
Michael Turner's avatar
Michael Turner committed
297
    {
298
        // we need to do the boundary layer generation in a face by face basis
Michael Turner's avatar
Michael Turner committed
299
300
        MakeBLPrep();

301
302
303
304
305
        // Im going to do a horrendous trick to get the edge orientaion.
        // Going to activate the first routine of facemeshing without actually
        // face meshing, this will orientate the edgeloop objects (hopefully);
        // which can be used by the makebl command to know the normal
        // orienation
Michael Turner's avatar
Michael Turner committed
306
307
308
        for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
        {
            m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
309
                i, m_mesh, m_curvemeshes, 100);
Michael Turner's avatar
Michael Turner committed
310
311
312
313

            m_facemeshes[i]->OrientateCurves();
            MakeBL(i, m_facemeshes[i]->GetEdges());
        }
314
315
    }

Julian Marcon's avatar
Tidy.    
Julian Marcon committed
316
    // m_mesh->m_element[1].clear();
317

318
319
320
321
322
323
324
325
326
327
328
329
330
    if (m_mesh->m_verbose)
    {
        cout << endl << "\tFace meshing:" << endl << endl;
    }

    // linear mesh all surfaces
    for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
    {
        if (m_mesh->m_verbose)
        {
            LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumSurf(),
                                           "Face progress");
        }
331

332
        m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
333
            i, m_mesh, m_curvemeshes, 100);
334

335
336
337
338
339
340
341
342
343
344
345
346
        m_facemeshes[i]->Mesh();
    }

    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
    ProcessComposites();

    Report();
}

Michael Turner's avatar
Michael Turner committed
347
void Generator2D::MakeBLPrep()
Michael Turner's avatar
Michael Turner committed
348
{
349
350
351
352
353
354
355
    if (m_mesh->m_verbose)
    {
        cout << endl << "\tBoundary layer meshing:" << endl << endl;
    }

    // identify the nodes which will become the boundary layer.
    ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
Michael Turner's avatar
Michael Turner committed
356
357
                                  m_blCurves);
    m_thickness = m_config["blthick"].as<NekDouble>();
358

359
360
    for (vector<unsigned>::iterator it = m_blCurves.begin();
         it != m_blCurves.end(); ++it)
361
    {
362
363
        vector<EdgeSharedPtr> localedges = m_curvemeshes[*it]->GetMeshEdges();
        for (int i = 0; i < localedges.size(); i++)
364
        {
Michael Turner's avatar
Michael Turner committed
365
366
            m_nodesToEdge[localedges[i]->m_n1].push_back(localedges[i]);
            m_nodesToEdge[localedges[i]->m_n2].push_back(localedges[i]);
367
368
        }
    }
Michael Turner's avatar
Michael Turner committed
369
}
370

Julian Marcon's avatar
Julian Marcon committed
371
void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
Michael Turner's avatar
Michael Turner committed
372
{
373
    map<int, int> edgeToOrient;
Julian Marcon's avatar
Julian Marcon committed
374
    for (vector<EdgeLoop>::iterator lit = e.begin(); lit != e.end(); ++lit)
375
    {
Julian Marcon's avatar
Julian Marcon committed
376
377
        for (int i = 0; i < lit->edges.size(); ++i)
        {
378
379
380
            edgeToOrient[lit->edges[i]->GetId()] = lit->edgeo[i];
        }
    }
381

382
    map<int, Array<OneD, NekDouble> > edgeNormals;
Julian Marcon's avatar
Julian Marcon committed
383

384
    int eid = 0;
Julian Marcon's avatar
Julian Marcon committed
385

386
    for (vector<unsigned>::iterator it = m_blCurves.begin();
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
387
         it != m_blCurves.end(); ++it)
388
389
    {
        int edgeo = edgeToOrient[*it];
Julian Marcon's avatar
Julian Marcon committed
390

391
        vector<EdgeSharedPtr> es = m_curvemeshes[*it]->GetMeshEdges();
Julian Marcon's avatar
Julian Marcon committed
392

393
394
        // on each !!!EDGE!!! calculate a normal
        // always to the left unless edgeo is 1
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
395
        for (int j = 0; j < es.size(); j++)
396
397
        {
            es[j]->m_id = eid++;
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
398
399
400
401
            Array<OneD, NekDouble> p1 =
                (edgeo == 0) ? es[j]->m_n1->GetLoc() : es[j]->m_n2->GetLoc();
            Array<OneD, NekDouble> p2 =
                (edgeo == 0) ? es[j]->m_n2->GetLoc() : es[j]->m_n1->GetLoc();
402
            Array<OneD, NekDouble> n(2);
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
403
404
405
            n[0]          = p1[1] - p2[1];
            n[1]          = p2[0] - p1[0];
            NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
406
407
408
409
410
            n[0] /= mag;
            n[1] /= mag;
            edgeNormals[es[j]->m_id] = n;
        }
    }
Julian Marcon's avatar
Julian Marcon committed
411

412
413
    map<NodeSharedPtr, NodeSharedPtr> nodeNormals;
    map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
414
    for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
415
416
    {
        Array<OneD, NekDouble> n(3);
Dave Moxey's avatar
Dave Moxey committed
417
418
        ASSERTL0(it->second.size() == 2,
                 "wierdness, most likely bl_surfs are incorrect");
419
420
421
        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
422
423
424
        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]);
425
426
427
428
429
430
431
        n[0] /= mag;
        n[1] /= mag;

        n[0] = n[0] * m_thickness + it->first->m_x;
        n[1] = n[1] * m_thickness + it->first->m_y;
        n[2] = 0.0;

Dave Moxey's avatar
Dave Moxey committed
432
433
        NodeSharedPtr nn = boost::shared_ptr<Node>(
            new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
434
435
        CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
        Array<OneD, NekDouble> uv = s->locuv(n);
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
436
        nn->SetCADSurf(faceid, s, uv);
437
438
        nodeNormals[it->first] = nn;
    }
439

440
    for (vector<unsigned>::iterator it = m_blCurves.begin();
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
441
         it != m_blCurves.end(); ++it)
442
443
    {
        int edgeo = edgeToOrient[*it];
Julian Marcon's avatar
Julian Marcon committed
444

445
446
        vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
        vector<NodeSharedPtr> newNs;
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
447
        for (int i = 0; i < ns.size(); i++)
448
449
450
        {
            newNs.push_back(nodeNormals[ns[i]]);
        }
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
451
452
        m_curvemeshes[*it] =
            MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);
Julian Marcon's avatar
Julian Marcon committed
453

Julian Marcon's avatar
Tidy.    
Julian Marcon committed
454
        if (edgeo == 1)
455
456
457
458
459
460
        {
            reverse(ns.begin(), ns.end());
        }
        for (int i = 0; i < ns.size() - 1; ++i)
        {
            vector<NodeSharedPtr> qns;
Julian Marcon's avatar
Julian Marcon committed
461

462
            qns.push_back(ns[i]);
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
463
464
            qns.push_back(ns[i + 1]);
            qns.push_back(nodeNormals[ns[i + 1]]);
465
            qns.push_back(nodeNormals[ns[i]]);
Julian Marcon's avatar
Julian Marcon committed
466

467
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
Julian Marcon's avatar
Julian Marcon committed
468

469
            vector<int> tags;
470
            tags.push_back(101);
Julian Marcon's avatar
Julian Marcon committed
471

472
473
            ElementSharedPtr E = GetElementFactory().CreateInstance(
                LibUtilities::eQuadrilateral, conf, qns, tags);
Julian Marcon's avatar
Julian Marcon committed
474

475
            E->m_parentCAD = m_mesh->m_cad->GetSurf(faceid);
476
477
478

            for (int j = 0; j < E->GetEdgeCount(); ++j)
            {
Julian Marcon's avatar
Tidy.    
Julian Marcon committed
479
                pair<EdgeSet::iterator, bool> testIns;
480
481
482
483
484
485
486
487
488
                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);
                }
            }
489
            m_mesh->m_element[2].push_back(E);
Julian Marcon's avatar
Julian Marcon committed
490
        }
491
    }
Michael Turner's avatar
Michael Turner committed
492
493
}

494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
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;
    }
}
509
510
}
}