CurveMesh.cpp 11.9 KB
Newer Older
Michael Turner's avatar
Michael Turner committed
1
2
////////////////////////////////////////////////////////////////////////////////
//
Michael Turner's avatar
Michael Turner committed
3
//  File: CurveMesh.cpp
Michael Turner's avatar
Michael Turner committed
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.
//
Michael Turner's avatar
Michael Turner committed
32
//  Description: curvemesh object methods.
Michael Turner's avatar
Michael Turner committed
33
34
35
//
////////////////////////////////////////////////////////////////////////////////

36
#include <NekMeshUtils/Octree/Octree.h>
Michael Turner's avatar
Michael Turner committed
37
#include <NekMeshUtils/SurfaceMeshing/CurveMesh.h>
Michael Turner's avatar
Michael Turner committed
38
39

using namespace std;
40
41
namespace Nektar
{
42
namespace NekMeshUtils
43
{
Michael Turner's avatar
Michael Turner committed
44

45
void CurveMesh::Mesh()
Michael Turner's avatar
Michael Turner committed
46
{
Michael Turner's avatar
Michael Turner committed
47
    // this algorithm is mostly based on the work in chapter 19
48

mike turner's avatar
mike turner committed
49
    m_bounds      = m_cadcurve->GetBounds();
Michael Turner's avatar
Michael Turner committed
50
    m_curvelength = m_cadcurve->GetTotLength();
Michael Turner's avatar
Michael Turner committed
51
52
53
    m_numSamplePoints =
        int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 10;
    ds = m_curvelength / (m_numSamplePoints - 1);
Michael Turner's avatar
Michael Turner committed
54

Julian Marcon's avatar
Julian Marcon committed
55
56
57
58
59
60
61
62
63
    NekDouble totalOffset = 0.0;
    for (map<unsigned, NekDouble>::iterator ie = m_endoffset.begin();
         ie != m_endoffset.end(); ++ie)
    {
        totalOffset += ie->second;
    }
    ASSERTL0(m_curvelength > totalOffset,
             "Boundary layers too thick for adjacent curve");

Michael Turner's avatar
Michael Turner committed
64
    GetSampleFunction();
Michael Turner's avatar
Michael Turner committed
65

Michael Turner's avatar
Michael Turner committed
66
    Ae = 0.0;
Michael Turner's avatar
Michael Turner committed
67

68
    for (int i = 0; i < m_numSamplePoints - 1; i++)
Michael Turner's avatar
Michael Turner committed
69
    {
70
        Ae += ds * (1.0 / m_dst[i][0] + 1.0 / m_dst[i + 1][0]) / 2.0;
Michael Turner's avatar
Michael Turner committed
71
    }
Michael Turner's avatar
Michael Turner committed
72

73
    Ne = round(Ae);
Michael Turner's avatar
Michael Turner committed
74

Julian Marcon's avatar
Julian Marcon committed
75
    if (Ne + 1 < 2 + m_endoffset.size())
Michael Turner's avatar
Michael Turner committed
76
    {
Julian Marcon's avatar
Julian Marcon committed
77
78
79
        Ne = 1 + m_endoffset.size();

        meshsvalue.resize(Ne + 1);
Michael Turner's avatar
Michael Turner committed
80
81
        meshsvalue[0] = 0.0;
        meshsvalue[1] = m_curvelength;
Julian Marcon's avatar
Julian Marcon committed
82
83
84
85
86
87
88
89
90

        if (m_endoffset.count(0))
        {
            meshsvalue[1] = m_endoffset[0];
        }
        if (m_endoffset.count(1))
        {
            meshsvalue[Ne - 1] = m_curvelength - m_endoffset[1];
        }
Michael Turner's avatar
Michael Turner committed
91
    }
Michael Turner's avatar
Michael Turner committed
92
    else
93
    {
Michael Turner's avatar
Michael Turner committed
94

Michael Turner's avatar
Michael Turner committed
95
        GetPhiFunction();
Michael Turner's avatar
Michael Turner committed
96

Michael Turner's avatar
Michael Turner committed
97
98
99
100
        meshsvalue.resize(Ne + 1);
        meshsvalue[0]  = 0.0;
        meshsvalue[Ne] = m_curvelength;

Julian Marcon's avatar
Julian Marcon committed
101
102
103
104
105
106
107
108
109
110
111
        if (m_endoffset.count(0))
        {
            meshsvalue[1] = m_endoffset[0];
        }
        if (m_endoffset.count(1))
        {
            meshsvalue[Ne - 1] = m_curvelength - m_endoffset[1];
        }

        for (int i = 1 + m_endoffset.count(0); i < Ne - m_endoffset.count(1);
             i++)
Michael Turner's avatar
Michael Turner committed
112
        {
113
114
115
            int iterationcounter = 0;
            bool iterate         = true;
            int k                = i;
Michael Turner's avatar
Michael Turner committed
116
            NekDouble ski        = meshsvalue[i - 1];
Michael Turner's avatar
Michael Turner committed
117
            NekDouble lastSki;
118
            while (iterate)
Michael Turner's avatar
Michael Turner committed
119
            {
Michael Turner's avatar
Michael Turner committed
120
                iterationcounter++;
121
122
                NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
                lastSki       = ski;
Michael Turner's avatar
Michael Turner committed
123
                ski           = ski - rhs;
124
                if (abs(lastSki - ski) < 1E-8)
Michael Turner's avatar
Michael Turner committed
125
                {
Michael Turner's avatar
Michael Turner committed
126
                    iterate = false;
Michael Turner's avatar
Michael Turner committed
127
                }
Michael Turner's avatar
Michael Turner committed
128

129
                ASSERTL0(iterationcounter < 1000000, "iteration failed");
Michael Turner's avatar
Michael Turner committed
130
            }
Michael Turner's avatar
Michael Turner committed
131

Michael Turner's avatar
Michael Turner committed
132
            meshsvalue[i] = ski;
Michael Turner's avatar
Michael Turner committed
133
        }
Michael Turner's avatar
Michael Turner committed
134
    }
Michael Turner's avatar
Michael Turner committed
135

Michael Turner's avatar
Michael Turner committed
136
137
    NekDouble t;
    Array<OneD, NekDouble> loc;
Michael Turner's avatar
Michael Turner committed
138

139
    vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
Julian Marcon's avatar
Julian Marcon committed
140
141
    vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s =
        m_cadcurve->GetAdjSurf();
Michael Turner's avatar
Michael Turner committed
142

143
    NodeSharedPtr n = verts[0]->GetNode();
Michael Turner's avatar
Michael Turner committed
144
    t               = m_bounds[0];
145
    n->SetCADCurve(m_id, m_cadcurve, t);
146
    loc = n->GetLoc();
147
    for (int j = 0; j < s.size(); j++)
148
    {
mike turner's avatar
mike turner committed
149
150
151
152
        if (verts[0]->IsDegen() == s[j].first->GetId())
        {
            // if the degen has been set for this node the node
            // already knows its corrected location
153
            continue;
mike turner's avatar
mike turner committed
154
        }
155

mike turner's avatar
mike turner committed
156
157
        Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
        n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
158
159
    }
    m_meshpoints.push_back(n);
Michael Turner's avatar
Michael Turner committed
160

161
    for (int i = 1; i < meshsvalue.size() - 1; i++)
Michael Turner's avatar
Michael Turner committed
162
    {
163
164
        t                = m_cadcurve->tAtArcLength(meshsvalue[i]);
        loc              = m_cadcurve->P(t);
165
        NodeSharedPtr n2 = boost::shared_ptr<Node>(
166
            new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
167
        n2->SetCADCurve(m_id, m_cadcurve, t);
168
        for (int j = 0; j < s.size(); j++)
169
        {
mike turner's avatar
mike turner committed
170
171
            Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
            n2->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
172
        }
173
        m_meshpoints.push_back(n2);
Michael Turner's avatar
Michael Turner committed
174
    }
Michael Turner's avatar
Michael Turner committed
175

176
    n = verts[1]->GetNode();
Michael Turner's avatar
Michael Turner committed
177
    t = m_bounds[1];
178
    n->SetCADCurve(m_id, m_cadcurve, t);
179
    loc = n->GetLoc();
180
    for (int j = 0; j < s.size(); j++)
181
    {
mike turner's avatar
mike turner committed
182
183
184
185
        if (verts[1]->IsDegen() == s[j].first->GetId())
        {
            // if the degen has been set for this node the node
            // already knows its corrected location
186
            continue;
mike turner's avatar
mike turner committed
187
        }
188

mike turner's avatar
mike turner committed
189
190
        Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
        n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
191
192
    }
    m_meshpoints.push_back(n);
Michael Turner's avatar
Michael Turner committed
193

194
195
    ASSERTL0(Ne + 1 == m_meshpoints.size(),
             "incorrect number of points in curve mesh");
196

Michael Turner's avatar
Michael Turner committed
197
198
199
200
201
202
203
204
205
    // make edges and add them to the edgeset for the face mesher to use
    for (int i = 0; i < m_meshpoints.size() - 1; i++)
    {
        EdgeSharedPtr e = boost::shared_ptr<Edge>(
            new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
        e->m_parentCAD = m_cadcurve;
        m_mesh->m_edgeSet.insert(e);
        m_meshedges.push_back(e);
    }
206

Michael Turner's avatar
Michael Turner committed
207
    if (m_mesh->m_verbose)
208
    {
Michael Turner's avatar
Michael Turner committed
209
210
        cout << "\r                                                            "
                "    "
211
212
213
214
215
216
217
                "                             ";
        cout << scientific << "\r\t\tCurve " << m_id << endl
             << "\t\t\tLength: " << m_curvelength << endl
             << "\t\t\tNodes: " << m_meshpoints.size() << endl
             << "\t\t\tSample points: " << m_numSamplePoints << endl
             << endl;
    }
218
219
}

Michael Turner's avatar
Michael Turner committed
220
221
222
223
224
void CurveMesh::GetPhiFunction()
{
    m_ps.resize(m_numSamplePoints);
    vector<NekDouble> newPhi;
    newPhi.resize(2);
Michael Turner's avatar
Michael Turner committed
225

226
    newPhi[0] = 0.0;
Michael Turner's avatar
Michael Turner committed
227
    newPhi[1] = 0.0;
Michael Turner's avatar
Michael Turner committed
228

229
    m_ps[0] = newPhi;
Michael Turner's avatar
Michael Turner committed
230

231
    NekDouble runningInt = 0.0;
Michael Turner's avatar
Michael Turner committed
232

233
    for (int i = 1; i < m_numSamplePoints; i++)
Michael Turner's avatar
Michael Turner committed
234
    {
235
236
237
238
        runningInt += (1.0 / m_dst[i - 1][0] + 1.0 / m_dst[i][0]) / 2.0 * ds;
        newPhi[0] = Ne / Ae * runningInt;
        newPhi[1] = m_dst[i][1];
        m_ps[i]   = newPhi;
Michael Turner's avatar
Michael Turner committed
239
    }
Michael Turner's avatar
Michael Turner committed
240
241
242
243
}

NekDouble CurveMesh::EvaluateDS(NekDouble s)
{
244
245
    int a = 0;
    int b = 0;
Michael Turner's avatar
Michael Turner committed
246

Julian Marcon's avatar
Julian Marcon committed
247
    ASSERTL1(!(s < 0) && !(s > m_curvelength), "s out of bounds");
Michael Turner's avatar
Michael Turner committed
248

Michael Turner's avatar
Michael Turner committed
249
    if (s == 0)
Michael Turner's avatar
Michael Turner committed
250
    {
Michael Turner's avatar
Michael Turner committed
251
252
        return m_dst[0][0];
    }
Michael Turner's avatar
Michael Turner committed
253
    else if (s == m_curvelength)
Michael Turner's avatar
Michael Turner committed
254
    {
255
        return m_dst[m_numSamplePoints - 1][0];
Michael Turner's avatar
Michael Turner committed
256
    }
Michael Turner's avatar
Michael Turner committed
257

258
    for (int i = 0; i < m_numSamplePoints - 1; i++)
Michael Turner's avatar
Michael Turner committed
259
    {
260
        if (m_dst[i][1] < s && m_dst[i + 1][1] >= s)
Michael Turner's avatar
Michael Turner committed
261
        {
262
263
            a = i;
            b = i + 1;
Michael Turner's avatar
Michael Turner committed
264
            break;
Michael Turner's avatar
Michael Turner committed
265
        }
Michael Turner's avatar
Michael Turner committed
266
    }
Michael Turner's avatar
Michael Turner committed
267

Michael Turner's avatar
Michael Turner committed
268
269
270
271
    NekDouble s1 = m_dst[a][1];
    NekDouble s2 = m_dst[b][1];
    NekDouble d1 = m_dst[a][0];
    NekDouble d2 = m_dst[b][0];
Michael Turner's avatar
Michael Turner committed
272

273
274
    NekDouble m = (d2 - d1) / (s2 - s1);
    NekDouble c = d2 - m * s2;
Michael Turner's avatar
Michael Turner committed
275

276
    ASSERTL0(m * s + c == m * s + c, "DS"); // was getting nans here
Michael Turner's avatar
Michael Turner committed
277

278
    return m * s + c;
Michael Turner's avatar
Michael Turner committed
279
}
Michael Turner's avatar
Michael Turner committed
280

Michael Turner's avatar
Michael Turner committed
281
282
NekDouble CurveMesh::EvaluatePS(NekDouble s)
{
283
284
    int a = 0;
    int b = 0;
Michael Turner's avatar
Michael Turner committed
285

Julian Marcon's avatar
Julian Marcon committed
286
    ASSERTL1(!(s < 0) && !(s > m_curvelength), "s out of bounds");
Michael Turner's avatar
Michael Turner committed
287

Michael Turner's avatar
Michael Turner committed
288
    if (s == 0)
Michael Turner's avatar
Michael Turner committed
289
    {
Michael Turner's avatar
Michael Turner committed
290
291
        return m_ps[0][0];
    }
Michael Turner's avatar
Michael Turner committed
292
    else if (s == m_curvelength)
Michael Turner's avatar
Michael Turner committed
293
    {
294
        return m_ps[m_numSamplePoints - 1][0];
Michael Turner's avatar
Michael Turner committed
295
    }
Michael Turner's avatar
Michael Turner committed
296

297
    for (int i = 0; i < m_numSamplePoints - 1; i++)
Michael Turner's avatar
Michael Turner committed
298
    {
299
        if (m_ps[i][1] < s && m_ps[i + 1][1] >= s)
Michael Turner's avatar
Michael Turner committed
300
        {
301
302
            a = i;
            b = i + 1;
Michael Turner's avatar
Michael Turner committed
303
            break;
Michael Turner's avatar
Michael Turner committed
304
        }
Michael Turner's avatar
Michael Turner committed
305
    }
Michael Turner's avatar
Michael Turner committed
306

307
    if (a == b)
Michael Turner's avatar
Michael Turner committed
308
309
310
311
312
313
    {
        cout << endl;
        cout << a << " " << b << endl;
        cout << s << endl;
        exit(-1);
    }
Michael Turner's avatar
Michael Turner committed
314

Michael Turner's avatar
Michael Turner committed
315
316
317
318
    NekDouble s1 = m_ps[a][1];
    NekDouble s2 = m_ps[b][1];
    NekDouble d1 = m_ps[a][0];
    NekDouble d2 = m_ps[b][0];
Michael Turner's avatar
Michael Turner committed
319

320
321
    NekDouble m = (d2 - d1) / (s2 - s1);
    NekDouble c = d2 - m * s2;
Michael Turner's avatar
Michael Turner committed
322

323
    ASSERTL0(m * s + c == m * s + c, "PS");
Michael Turner's avatar
Michael Turner committed
324

325
    return m * s + c;
Michael Turner's avatar
Michael Turner committed
326
}
Michael Turner's avatar
Michael Turner committed
327

Michael Turner's avatar
Michael Turner committed
328
329
330
void CurveMesh::GetSampleFunction()
{
    m_dst.resize(m_numSamplePoints);
Michael Turner's avatar
Michael Turner committed
331

Michael Turner's avatar
Michael Turner committed
332
333
    vector<NekDouble> dsti;
    dsti.resize(3);
Michael Turner's avatar
Michael Turner committed
334

335
    for (int i = 0; i < m_numSamplePoints; i++)
Michael Turner's avatar
Michael Turner committed
336
    {
Michael Turner's avatar
Michael Turner committed
337
338
        dsti[1]     = i * ds;
        NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
Michael Turner's avatar
Michael Turner committed
339

Michael Turner's avatar
Michael Turner committed
340
        Array<OneD, NekDouble> loc = m_cadcurve->P(t);
Michael Turner's avatar
Michael Turner committed
341

Julian Marcon's avatar
Julian Marcon committed
342
        bool found = false;
Michael Turner's avatar
Michael Turner committed
343

Julian Marcon's avatar
Julian Marcon committed
344
        if (m_endoffset.count(0))
345
        {
Julian Marcon's avatar
Julian Marcon committed
346
            if (dsti[1] < m_endoffset[0])
Michael Turner's avatar
Michael Turner committed
347
            {
Julian Marcon's avatar
Julian Marcon committed
348
                dsti[0] = m_endoffset[0];
Julian Marcon's avatar
Julian Marcon committed
349
                found   = true;
Michael Turner's avatar
Michael Turner committed
350
            }
Julian Marcon's avatar
Julian Marcon committed
351
352
353
354
        }
        if (m_endoffset.count(1) && !found)
        {
            if (dsti[1] > m_curvelength - m_endoffset[1])
Michael Turner's avatar
Michael Turner committed
355
            {
Julian Marcon's avatar
Julian Marcon committed
356
                dsti[0] = m_endoffset[1];
Julian Marcon's avatar
Julian Marcon committed
357
                found   = true;
Michael Turner's avatar
Michael Turner committed
358
            }
359
        }
Julian Marcon's avatar
Julian Marcon committed
360
361
        if (!found)
        {
Michael Turner's avatar
Michael Turner committed
362
            dsti[0] = m_mesh->m_octree->Query(loc);
Julian Marcon's avatar
Julian Marcon committed
363
        }
364

365
        dsti[2] = t;
Michael Turner's avatar
Michael Turner committed
366

367
        m_dst[i] = dsti;
Michael Turner's avatar
Michael Turner committed
368
    }
Michael Turner's avatar
Michael Turner committed
369
}
Michael Turner's avatar
Michael Turner committed
370
371
372

void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
{
Julian Marcon's avatar
Julian Marcon committed
373
    // clear current mesh points and remove edges from edgeset
Michael Turner's avatar
Michael Turner committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
    m_meshpoints.clear();
    for (int i = 0; i < m_meshedges.size(); i++)
    {
        m_mesh->m_edgeSet.erase(m_meshedges[i]);
    }
    m_meshedges.clear();

    ///////

    int tid = from->GetId();
    Array<OneD, NekDouble> T =
        m_mesh->m_cad->GetPeriodicTranslationVector(tid, m_id);

    CADCurveSharedPtr c1 = m_mesh->m_cad->GetCurve(tid);

    bool reversed = c1->GetOrienationWRT(1) == m_cadcurve->GetOrienationWRT(1);

    vector<NodeSharedPtr> nodes = from->GetMeshPoints();

    vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > surfs =
        m_cadcurve->GetAdjSurf();

    for (int i = 1; i < nodes.size() - 1; i++)
    {
        Array<OneD, NekDouble> loc = nodes[i]->GetLoc();
Julian Marcon's avatar
Julian Marcon committed
399
400
        NodeSharedPtr nn = NodeSharedPtr(
            new Node(m_mesh->m_numNodes++, loc[0] + T[0], loc[1] + T[1], 0.0));
Michael Turner's avatar
Michael Turner committed
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422

        for (int j = 0; j < surfs.size(); j++)
        {
            nn->SetCADSurf(surfs[j].first->GetId(), surfs[j].first,
                           surfs[j].first->locuv(nn->GetLoc()));
        }

        nn->SetCADCurve(m_id, m_cadcurve, m_cadcurve->loct(nn->GetLoc()));

        m_meshpoints.push_back(nn);
    }

    // Reverse internal nodes of the vector if necessary
    if (reversed)
    {
        reverse(m_meshpoints.begin(), m_meshpoints.end());
    }

    vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();

    m_meshpoints.insert(m_meshpoints.begin(), verts[0]->GetNode());
    m_meshpoints.push_back(verts[1]->GetNode());
Julian Marcon's avatar
Julian Marcon committed
423
    // dont need to realign cad for vertices
Michael Turner's avatar
Michael Turner committed
424
425
426
427
428
429
430
431
432
433
434

    // make edges and add them to the edgeset for the face mesher to use
    for (int i = 0; i < m_meshpoints.size() - 1; i++)
    {
        EdgeSharedPtr e = boost::shared_ptr<Edge>(
            new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
        e->m_parentCAD = m_cadcurve;
        m_mesh->m_edgeSet.insert(e);
        m_meshedges.push_back(e);
    }
}
Michael Turner's avatar
Michael Turner committed
435
436
}
}