CurveMesh.cpp 11.4 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

Michael Turner's avatar
Michael Turner committed
75
    if (Ne + 1 < 2)
Michael Turner's avatar
Michael Turner committed
76
    {
Michael Turner's avatar
Michael Turner committed
77
78
79
80
        meshsvalue.resize(2);
        meshsvalue[0] = 0.0;
        meshsvalue[1] = m_curvelength;
        Ne            = 1;
Michael Turner's avatar
Michael Turner committed
81
    }
Michael Turner's avatar
Michael Turner committed
82
    else
83
    {
Michael Turner's avatar
Michael Turner committed
84

Michael Turner's avatar
Michael Turner committed
85
        GetPhiFunction();
Michael Turner's avatar
Michael Turner committed
86

Michael Turner's avatar
Michael Turner committed
87
88
89
90
        meshsvalue.resize(Ne + 1);
        meshsvalue[0]  = 0.0;
        meshsvalue[Ne] = m_curvelength;

91
        for (int i = 1; i < Ne; i++)
Michael Turner's avatar
Michael Turner committed
92
        {
93
94
95
            int iterationcounter = 0;
            bool iterate         = true;
            int k                = i;
Michael Turner's avatar
Michael Turner committed
96
            NekDouble ski        = meshsvalue[i - 1];
Michael Turner's avatar
Michael Turner committed
97
            NekDouble lastSki;
98
            while (iterate)
Michael Turner's avatar
Michael Turner committed
99
            {
Michael Turner's avatar
Michael Turner committed
100
                iterationcounter++;
101
102
                NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
                lastSki       = ski;
Michael Turner's avatar
Michael Turner committed
103
                ski           = ski - rhs;
104
                if (abs(lastSki - ski) < 1E-8)
Michael Turner's avatar
Michael Turner committed
105
                {
Michael Turner's avatar
Michael Turner committed
106
                    iterate = false;
Michael Turner's avatar
Michael Turner committed
107
                }
Michael Turner's avatar
Michael Turner committed
108

109
                ASSERTL0(iterationcounter < 1000000, "iteration failed");
Michael Turner's avatar
Michael Turner committed
110
            }
Michael Turner's avatar
Michael Turner committed
111

Michael Turner's avatar
Michael Turner committed
112
            meshsvalue[i] = ski;
Michael Turner's avatar
Michael Turner committed
113
        }
Michael Turner's avatar
Michael Turner committed
114
    }
Michael Turner's avatar
Michael Turner committed
115

Michael Turner's avatar
Michael Turner committed
116
117
    NekDouble t;
    Array<OneD, NekDouble> loc;
Michael Turner's avatar
Michael Turner committed
118

119
    vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
Julian Marcon's avatar
Julian Marcon committed
120
121
    vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s =
        m_cadcurve->GetAdjSurf();
Michael Turner's avatar
Michael Turner committed
122

123
    NodeSharedPtr n = verts[0]->GetNode();
Michael Turner's avatar
Michael Turner committed
124
    t               = m_bounds[0];
125
    n->SetCADCurve(m_id, m_cadcurve, t);
126
    loc = n->GetLoc();
127
    for (int j = 0; j < s.size(); j++)
128
    {
mike turner's avatar
mike turner committed
129
130
131
132
        if (verts[0]->IsDegen() == s[j].first->GetId())
        {
            // if the degen has been set for this node the node
            // already knows its corrected location
133
            continue;
mike turner's avatar
mike turner committed
134
        }
135

mike turner's avatar
mike turner committed
136
137
        Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
        n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
138
139
    }
    m_meshpoints.push_back(n);
Michael Turner's avatar
Michael Turner committed
140

141
    for (int i = 1; i < meshsvalue.size() - 1; i++)
Michael Turner's avatar
Michael Turner committed
142
    {
143
144
        t                = m_cadcurve->tAtArcLength(meshsvalue[i]);
        loc              = m_cadcurve->P(t);
145
        NodeSharedPtr n2 = boost::shared_ptr<Node>(
146
            new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
147
        n2->SetCADCurve(m_id, m_cadcurve, t);
148
        for (int j = 0; j < s.size(); j++)
149
        {
mike turner's avatar
mike turner committed
150
151
            Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
            n2->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
152
        }
153
        m_meshpoints.push_back(n2);
Michael Turner's avatar
Michael Turner committed
154
    }
Michael Turner's avatar
Michael Turner committed
155

156
    n = verts[1]->GetNode();
Michael Turner's avatar
Michael Turner committed
157
    t = m_bounds[1];
158
    n->SetCADCurve(m_id, m_cadcurve, t);
159
    loc = n->GetLoc();
160
    for (int j = 0; j < s.size(); j++)
161
    {
mike turner's avatar
mike turner committed
162
163
164
165
        if (verts[1]->IsDegen() == s[j].first->GetId())
        {
            // if the degen has been set for this node the node
            // already knows its corrected location
166
            continue;
mike turner's avatar
mike turner committed
167
        }
168

mike turner's avatar
mike turner committed
169
170
        Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
        n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
171
172
    }
    m_meshpoints.push_back(n);
Michael Turner's avatar
Michael Turner committed
173

174
175
    ASSERTL0(Ne + 1 == m_meshpoints.size(),
             "incorrect number of points in curve mesh");
176

Michael Turner's avatar
Michael Turner committed
177
178
179
180
181
182
183
184
185
    // 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);
    }
186

Michael Turner's avatar
Michael Turner committed
187
    if (m_mesh->m_verbose)
188
    {
Michael Turner's avatar
Michael Turner committed
189
190
        cout << "\r                                                            "
                "    "
191
192
193
194
195
196
197
                "                             ";
        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;
    }
198
199
}

Michael Turner's avatar
Michael Turner committed
200
201
202
203
204
void CurveMesh::GetPhiFunction()
{
    m_ps.resize(m_numSamplePoints);
    vector<NekDouble> newPhi;
    newPhi.resize(2);
Michael Turner's avatar
Michael Turner committed
205

206
    newPhi[0] = 0.0;
Michael Turner's avatar
Michael Turner committed
207
    newPhi[1] = 0.0;
Michael Turner's avatar
Michael Turner committed
208

209
    m_ps[0] = newPhi;
Michael Turner's avatar
Michael Turner committed
210

211
    NekDouble runningInt = 0.0;
Michael Turner's avatar
Michael Turner committed
212

213
    for (int i = 1; i < m_numSamplePoints; i++)
Michael Turner's avatar
Michael Turner committed
214
    {
215
216
217
218
        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
219
    }
Michael Turner's avatar
Michael Turner committed
220
221
222
223
}

NekDouble CurveMesh::EvaluateDS(NekDouble s)
{
224
225
    int a = 0;
    int b = 0;
Michael Turner's avatar
Michael Turner committed
226

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

Michael Turner's avatar
Michael Turner committed
229
    if (s == 0)
Michael Turner's avatar
Michael Turner committed
230
    {
Michael Turner's avatar
Michael Turner committed
231
232
        return m_dst[0][0];
    }
Michael Turner's avatar
Michael Turner committed
233
    else if (s == m_curvelength)
Michael Turner's avatar
Michael Turner committed
234
    {
235
        return m_dst[m_numSamplePoints - 1][0];
Michael Turner's avatar
Michael Turner committed
236
    }
Michael Turner's avatar
Michael Turner committed
237

238
    for (int i = 0; i < m_numSamplePoints - 1; i++)
Michael Turner's avatar
Michael Turner committed
239
    {
240
        if (m_dst[i][1] < s && m_dst[i + 1][1] >= s)
Michael Turner's avatar
Michael Turner committed
241
        {
242
243
            a = i;
            b = i + 1;
Michael Turner's avatar
Michael Turner committed
244
            break;
Michael Turner's avatar
Michael Turner committed
245
        }
Michael Turner's avatar
Michael Turner committed
246
    }
Michael Turner's avatar
Michael Turner committed
247

Michael Turner's avatar
Michael Turner committed
248
249
250
251
    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
252

253
254
    NekDouble m = (d2 - d1) / (s2 - s1);
    NekDouble c = d2 - m * s2;
Michael Turner's avatar
Michael Turner committed
255

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

258
    return m * s + c;
Michael Turner's avatar
Michael Turner committed
259
}
Michael Turner's avatar
Michael Turner committed
260

Michael Turner's avatar
Michael Turner committed
261
262
NekDouble CurveMesh::EvaluatePS(NekDouble s)
{
263
264
    int a = 0;
    int b = 0;
Michael Turner's avatar
Michael Turner committed
265

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

Michael Turner's avatar
Michael Turner committed
268
    if (s == 0)
Michael Turner's avatar
Michael Turner committed
269
    {
Michael Turner's avatar
Michael Turner committed
270
271
        return m_ps[0][0];
    }
Michael Turner's avatar
Michael Turner committed
272
    else if (s == m_curvelength)
Michael Turner's avatar
Michael Turner committed
273
    {
274
        return m_ps[m_numSamplePoints - 1][0];
Michael Turner's avatar
Michael Turner committed
275
    }
Michael Turner's avatar
Michael Turner committed
276

277
    for (int i = 0; i < m_numSamplePoints - 1; i++)
Michael Turner's avatar
Michael Turner committed
278
    {
279
        if (m_ps[i][1] < s && m_ps[i + 1][1] >= s)
Michael Turner's avatar
Michael Turner committed
280
        {
281
282
            a = i;
            b = i + 1;
Michael Turner's avatar
Michael Turner committed
283
            break;
Michael Turner's avatar
Michael Turner committed
284
        }
Michael Turner's avatar
Michael Turner committed
285
    }
Michael Turner's avatar
Michael Turner committed
286

287
    if (a == b)
Michael Turner's avatar
Michael Turner committed
288
289
290
291
292
293
    {
        cout << endl;
        cout << a << " " << b << endl;
        cout << s << endl;
        exit(-1);
    }
Michael Turner's avatar
Michael Turner committed
294

Michael Turner's avatar
Michael Turner committed
295
296
297
298
    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
299

300
301
    NekDouble m = (d2 - d1) / (s2 - s1);
    NekDouble c = d2 - m * s2;
Michael Turner's avatar
Michael Turner committed
302

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

305
    return m * s + c;
Michael Turner's avatar
Michael Turner committed
306
}
Michael Turner's avatar
Michael Turner committed
307

Michael Turner's avatar
Michael Turner committed
308
309
310
void CurveMesh::GetSampleFunction()
{
    m_dst.resize(m_numSamplePoints);
Michael Turner's avatar
Michael Turner committed
311

Michael Turner's avatar
Michael Turner committed
312
313
    vector<NekDouble> dsti;
    dsti.resize(3);
Michael Turner's avatar
Michael Turner committed
314

315
    for (int i = 0; i < m_numSamplePoints; i++)
Michael Turner's avatar
Michael Turner committed
316
    {
Michael Turner's avatar
Michael Turner committed
317
318
        dsti[1]     = i * ds;
        NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
Michael Turner's avatar
Michael Turner committed
319

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

Julian Marcon's avatar
Julian Marcon committed
322
        bool found = false;
Michael Turner's avatar
Michael Turner committed
323

Julian Marcon's avatar
Julian Marcon committed
324
        if (m_endoffset.count(0))
325
        {
Julian Marcon's avatar
Julian Marcon committed
326
            if (dsti[1] < m_endoffset[0])
Michael Turner's avatar
Michael Turner committed
327
            {
Julian Marcon's avatar
Julian Marcon committed
328
329
                dsti[0] = m_endoffset[0];
                found = true;
Michael Turner's avatar
Michael Turner committed
330
            }
Julian Marcon's avatar
Julian Marcon committed
331
332
333
334
        }
        if (m_endoffset.count(1) && !found)
        {
            if (dsti[1] > m_curvelength - m_endoffset[1])
Michael Turner's avatar
Michael Turner committed
335
            {
Julian Marcon's avatar
Julian Marcon committed
336
337
                dsti[0] = m_endoffset[1];
                found = true;
Michael Turner's avatar
Michael Turner committed
338
            }
339
        }
Julian Marcon's avatar
Julian Marcon committed
340
341
        if (!found)
        {
Michael Turner's avatar
Michael Turner committed
342
            dsti[0] = m_mesh->m_octree->Query(loc);
Julian Marcon's avatar
Julian Marcon committed
343
        }
344

345
        dsti[2] = t;
Michael Turner's avatar
Michael Turner committed
346

347
        m_dst[i] = dsti;
Michael Turner's avatar
Michael Turner committed
348
    }
Michael Turner's avatar
Michael Turner committed
349
}
Michael Turner's avatar
Michael Turner committed
350
351
352

void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
{
Julian Marcon's avatar
Julian Marcon committed
353
    // clear current mesh points and remove edges from edgeset
Michael Turner's avatar
Michael Turner committed
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    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
379
380
        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
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402

        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
403
    // dont need to realign cad for vertices
Michael Turner's avatar
Michael Turner committed
404
405
406
407
408
409
410
411
412
413
414

    // 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
415
416
}
}