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