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

36
#include <NekMeshUtils/CADSystem/OCE/CADSurfOCE.h>
37
38

using namespace std;
39

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

45
46
std::string CADSurfOCE::key = GetCADSurfFactory().RegisterCreatorFunction(
    "oce", CADSurfOCE::create, "CADSurfOCE");
mike's avatar
mike committed
47

mike turner's avatar
mike turner committed
48
void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
Michael Turner's avatar
Michael Turner committed
49
{
50
51
    // this bit of code changes the units of the cad from mm opencascade
    // defualt to m
52
53

    m_s = BRep_Tool::Surface(TopoDS::Face(in));
Michael Turner's avatar
Michael Turner committed
54

Michael Turner's avatar
Michael Turner committed
55
    if (in.Orientation() == 1)
56
    {
57
        m_orientation = CADOrientation::eBackwards;
Michael Turner's avatar
Michael Turner committed
58
    }
59

Michael Turner's avatar
Michael Turner committed
60
61
62
63
64
    gp_Trsf transform;
    gp_Pnt ori(0.0, 0.0, 0.0);
    transform.SetScale(ori, 1.0 / 1000.0);
    TopLoc_Location mv(transform);
    in.Move(mv);
Michael Turner's avatar
Michael Turner committed
65

Michael Turner's avatar
Michael Turner committed
66
67
68
69
70
71
72
73
    m_occSurface = BRepAdaptor_Surface(TopoDS::Face(in));
    m_id         = i;

    m_bounds = Array<OneD, NekDouble>(4);
    BRepTools::UVBounds(TopoDS::Face(in), m_bounds[0], m_bounds[1], m_bounds[2],
                        m_bounds[3]);
    m_sas = new ShapeAnalysis_Surface(m_s);
    m_sas->SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
Michael Turner's avatar
Michael Turner committed
74
}
75

76
77
Array<OneD, NekDouble> CADSurfOCE::GetBounds()
{
Michael Turner's avatar
Michael Turner committed
78
    return m_bounds;
79
80
}

mike's avatar
mike committed
81
Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
Michael Turner's avatar
Michael Turner committed
82
{
83
    // has to transfer back to mm
Michael Turner's avatar
Michael Turner committed
84
    gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
85

86
    Array<OneD, NekDouble> uvr(2);
87

Michael Turner's avatar
Michael Turner committed
88
    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
Michael Turner's avatar
Michael Turner committed
89
90
    uvr[0]      = p2.X();
    uvr[1]      = p2.Y();
91

Michael Turner's avatar
Michael Turner committed
92
    gp_Pnt p3 = m_sas->Value(p2);
Michael Turner's avatar
Michael Turner committed
93
    if (p3.Distance(loc) > 1e-6)
Michael Turner's avatar
Michael Turner committed
94
    {
Michael Turner's avatar
Michael Turner committed
95
        cout << "large locuv distance " << p3.Distance(loc)/1000.0 << " " << m_id
Michael Turner's avatar
Michael Turner committed
96
             << endl;
97
    }
98

99
    // if the uv returned is slightly off the surface
Michael Turner's avatar
Michael Turner committed
100
    //(which ShapeAnalysis_Surface can do sometimes)
Michael Turner's avatar
Michael Turner committed
101
102
    if (uvr[0] < m_bounds[0] || uvr[0] > m_bounds[1] || uvr[1] < m_bounds[2] ||
        uvr[1] > m_bounds[3])
103
    {
Michael Turner's avatar
Michael Turner committed
104
        if (uvr[0] < m_bounds[0])
105
        {
Michael Turner's avatar
Michael Turner committed
106
            uvr[0] = m_bounds[0];
107
        }
Michael Turner's avatar
Michael Turner committed
108
        else if (uvr[0] > m_bounds[1])
109
        {
Michael Turner's avatar
Michael Turner committed
110
            uvr[0] = m_bounds[1];
111
        }
Michael Turner's avatar
Michael Turner committed
112
        else if (uvr[1] < m_bounds[2])
113
        {
Michael Turner's avatar
Michael Turner committed
114
            uvr[1] = m_bounds[2];
115
        }
Michael Turner's avatar
Michael Turner committed
116
        else if (uvr[1] > m_bounds[3])
117
        {
Michael Turner's avatar
Michael Turner committed
118
            uvr[1] = m_bounds[3];
119
120
121
        }
        else
        {
122
            ASSERTL0(false, "Cannot correct locuv");
123
        }
124
    }
125

Michael Turner's avatar
Michael Turner committed
126
127
    return uvr;
}
128

mike's avatar
mike committed
129
NekDouble CADSurfOCE::Curvature(Array<OneD, NekDouble> uv)
130
{
Michael Turner's avatar
Michael Turner committed
131
#if defined(NEKTAR_DEBUG)
132
    Test(uv);
Michael Turner's avatar
Michael Turner committed
133
#endif
134
135
136

    Array<OneD, NekDouble> n = N(uv);

137
138
139
    // a zero normal occurs at a signularity, CurvaturePoint
    // cannot be sampled here
    if (n[0] == 0 && n[1] == 0 && n[2] == 0)
140
141
142
143
144
145
    {
        return 0.0;
    }

    Array<OneD, NekDouble> r = D2(uv);

146
147
148
149
150
151
152
    // metric and curvature tensors
    NekDouble E = r[3] * r[3] + r[4] * r[4] + r[5] * r[5];
    NekDouble F = r[3] * r[6] + r[4] * r[7] + r[5] * r[8];
    NekDouble G = r[6] * r[6] + r[7] * r[7] + r[8] * r[8];
    NekDouble e = n[0] * r[9] + n[1] * r[10] + n[2] * r[11];
    NekDouble f = n[0] * r[15] + n[1] * r[16] + n[2] * r[17];
    NekDouble g = n[0] * r[12] + n[1] * r[13] + n[2] * r[14];
153

154
155
    // if det is zero cannot invert matrix, R=0 so must skip
    if (E * G - F * F < 1E-30)
156
    {
Michael Turner's avatar
Michael Turner committed
157
        return 0.0;
158
159
160
161
    }

    NekDouble K, H;

162
163
    K = (e * g - f * f) / (E * G - F * F);
    H = 0.5 * (e * G - 2 * f * F + g * E) / (E * G - F * F);
164
165

    NekDouble kv[2];
166
167
    kv[0] = abs(H + sqrt(H * H - K));
    kv[1] = abs(H - sqrt(H * H - K));
168
169
170
171

    return kv[0] > kv[1] ? kv[0] : kv[1];
}

mike's avatar
mike committed
172
NekDouble CADSurfOCE::DistanceTo(Array<OneD, NekDouble> p)
173
174
175
{
    gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);

Michael Turner's avatar
Michael Turner committed
176
    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
177

Michael Turner's avatar
Michael Turner committed
178
    gp_Pnt p3 = m_sas->Value(p2);
179
180
181
182

    return p3.Distance(loc);
}

Michael Turner's avatar
Michael Turner committed
183
184
void CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
                           Array<OneD, NekDouble> &uv)
185
186
187
{
    gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);

Michael Turner's avatar
Michael Turner committed
188
    gp_Pnt2d p2 = m_sas->ValueOfUV(loc, Precision::Confusion());
189

Michael Turner's avatar
Michael Turner committed
190
    gp_Pnt p3 = m_sas->Value(p2);
191

192
193
194
    tp[0] = p3.X() / 1000.0;
    tp[1] = p3.Y() / 1000.0;
    tp[2] = p3.Z() / 1000.0;
195
196

    uv[0] = p2.X();
Michael Turner's avatar
Michael Turner committed
197
    uv[1] = p2.Y();
198
199
}

mike's avatar
mike committed
200
Array<OneD, NekDouble> CADSurfOCE::P(Array<OneD, NekDouble> uv)
Michael Turner's avatar
Michael Turner committed
201
{
Michael Turner's avatar
Michael Turner committed
202
#if defined(NEKTAR_DEBUG)
Michael Turner's avatar
tweaks    
Michael Turner committed
203
    Test(uv);
Michael Turner's avatar
Michael Turner committed
204
#endif
Michael Turner's avatar
Michael Turner committed
205

Michael Turner's avatar
Michael Turner committed
206
207
    Array<OneD, NekDouble> location(3);
    gp_Pnt loc;
208
    loc         = m_occSurface.Value(uv[0], uv[1]);
Michael Turner's avatar
Michael Turner committed
209
210
211
212
213
    location[0] = loc.X();
    location[1] = loc.Y();
    location[2] = loc.Z();
    return location;
}
Michael Turner's avatar
Michael Turner committed
214

mike's avatar
mike committed
215
Array<OneD, NekDouble> CADSurfOCE::N(Array<OneD, NekDouble> uv)
Michael Turner's avatar
Michael Turner committed
216
{
Michael Turner's avatar
Michael Turner committed
217
#if defined(NEKTAR_DEBUG)
Michael Turner's avatar
tweaks    
Michael Turner committed
218
    Test(uv);
Michael Turner's avatar
Michael Turner committed
219
#endif
Michael Turner's avatar
Michael Turner committed
220

Michael Turner's avatar
Michael Turner committed
221
    BRepLProp_SLProps slp(m_occSurface, 2, 1e-8);
Michael Turner's avatar
Michael Turner committed
222
    slp.SetParameters(uv[0], uv[1]);
223

Michael Turner's avatar
Michael Turner committed
224
    if (!slp.IsNormalDefined())
Michael Turner's avatar
changes    
Michael Turner committed
225
    {
Michael Turner's avatar
Michael Turner committed
226
        return Array<OneD, NekDouble>(3, 0.0);
Michael Turner's avatar
changes    
Michael Turner committed
227
228
    }

Michael Turner's avatar
Michael Turner committed
229
230
231
232
    gp_Dir d = slp.Normal();

    Array<OneD, NekDouble> normal(3);

233
    if (m_orientation == CADOrientation::eBackwards)
Michael Turner's avatar
Michael Turner committed
234
    {
Michael Turner's avatar
Michael Turner committed
235
        d.Reverse();
236
    }
Michael Turner's avatar
Michael Turner committed
237

Michael Turner's avatar
Michael Turner committed
238
239
240
241
    normal[0] = d.X();
    normal[1] = d.Y();
    normal[2] = d.Z();

Michael Turner's avatar
Michael Turner committed
242
243
    return normal;
}
Michael Turner's avatar
Michael Turner committed
244

mike's avatar
mike committed
245
Array<OneD, NekDouble> CADSurfOCE::D1(Array<OneD, NekDouble> uv)
Michael Turner's avatar
Michael Turner committed
246
{
Michael Turner's avatar
Michael Turner committed
247
#if defined(NEKTAR_DEBUG)
Michael Turner's avatar
tweaks    
Michael Turner committed
248
    Test(uv);
Michael Turner's avatar
Michael Turner committed
249
#endif
Michael Turner's avatar
Michael Turner committed
250

Michael Turner's avatar
Michael Turner committed
251
252
253
    Array<OneD, NekDouble> r(9);
    gp_Pnt Loc;
    gp_Vec D1U, D1V;
254
    m_occSurface.D1(uv[0], uv[1], Loc, D1U, D1V);
Michael Turner's avatar
Michael Turner committed
255

256
257
258
259
260
261
262
263
264
    r[0] = Loc.X(); // x
    r[1] = Loc.Y(); // y
    r[2] = Loc.Z(); // z
    r[3] = D1U.X(); // dx/du
    r[4] = D1U.Y(); // dy/du
    r[5] = D1U.Z(); // dz/du
    r[6] = D1V.X(); // dx/dv
    r[7] = D1V.Y(); // dy/dv
    r[8] = D1V.Z(); // dz/dv
Michael Turner's avatar
Michael Turner committed
265
266
267

    return r;
}
Michael Turner's avatar
Michael Turner committed
268

mike's avatar
mike committed
269
Array<OneD, NekDouble> CADSurfOCE::D2(Array<OneD, NekDouble> uv)
Michael Turner's avatar
Michael Turner committed
270
{
Michael Turner's avatar
Michael Turner committed
271
#if defined(NEKTAR_DEBUG)
Michael Turner's avatar
tweaks    
Michael Turner committed
272
    Test(uv);
Michael Turner's avatar
Michael Turner committed
273
#endif
Michael Turner's avatar
Michael Turner committed
274

Michael Turner's avatar
Michael Turner committed
275
276
277
    Array<OneD, NekDouble> r(18);
    gp_Pnt Loc;
    gp_Vec D1U, D1V, D2U, D2V, D2UV;
278
    m_occSurface.D2(uv[0], uv[1], Loc, D1U, D1V, D2U, D2V, D2UV);
Michael Turner's avatar
Michael Turner committed
279

Michael Turner's avatar
Michael Turner committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    r[0]  = Loc.X();  // x
    r[1]  = Loc.Y();  // y
    r[2]  = Loc.Z();  // z
    r[3]  = D1U.X();  // dx/dx
    r[4]  = D1U.Y();  // dy/dy
    r[5]  = D1U.Z();  // dz/dz
    r[6]  = D1V.X();  // dx/dx
    r[7]  = D1V.Y();  // dy/dy
    r[8]  = D1V.Z();  // dz/dz
    r[9]  = D2U.X();  // d2x/du2
    r[10] = D2U.Y();  // d2y/du2
    r[11] = D2U.Z();  // d2z/du2
    r[12] = D2V.X();  // d2x/dv2
    r[13] = D2V.Y();  // d2y/dv2
    r[14] = D2V.Z();  // d2z/dv2
295
296
297
    r[15] = D2UV.X(); // d2x/dudv
    r[16] = D2UV.Y(); // d2y/dudv
    r[17] = D2UV.Z(); // d2z/dudv
Michael Turner's avatar
Michael Turner committed
298
299
300

    return r;
}
Michael Turner's avatar
Michael Turner committed
301

mike's avatar
mike committed
302
void CADSurfOCE::Test(Array<OneD, NekDouble> uv)
Michael Turner's avatar
tweaks    
Michael Turner committed
303
{
Dave Moxey's avatar
Dave Moxey committed
304
    stringstream error;
Michael Turner's avatar
tweaks    
Michael Turner committed
305
306
307
308
309

    error << "Point not within parameter plane: ";

    bool passed = true;

Michael Turner's avatar
Michael Turner committed
310
    if (uv[0] < m_bounds[0])
Michael Turner's avatar
tweaks    
Michael Turner committed
311
    {
Michael Turner's avatar
Michael Turner committed
312
        if (fabs(uv[0] - m_bounds[0]) > 1E-6)
313
        {
Michael Turner's avatar
Michael Turner committed
314
315
            error << "U(" << uv[0] << ") is less than Umin(" << m_bounds[0]
                  << ")";
316
317
            passed = false;
        }
Michael Turner's avatar
tweaks    
Michael Turner committed
318
    }
Michael Turner's avatar
Michael Turner committed
319
    else if (uv[0] > m_bounds[1])
Michael Turner's avatar
tweaks    
Michael Turner committed
320
    {
Michael Turner's avatar
Michael Turner committed
321
        if (fabs(uv[0] - m_bounds[1]) > 1E-6)
322
        {
Michael Turner's avatar
Michael Turner committed
323
324
            error << "U(" << uv[0] << ") is greater than Umax(" << m_bounds[1]
                  << ")";
325
326
            passed = false;
        }
Michael Turner's avatar
tweaks    
Michael Turner committed
327
    }
Michael Turner's avatar
Michael Turner committed
328
    else if (uv[1] < m_bounds[2])
Michael Turner's avatar
tweaks    
Michael Turner committed
329
    {
Michael Turner's avatar
Michael Turner committed
330
        if (fabs(uv[1] - m_bounds[2]) > 1E-6)
331
        {
Michael Turner's avatar
Michael Turner committed
332
333
            error << "V(" << uv[1] << ") is less than Vmin(" << m_bounds[2]
                  << ")";
334
335
            passed = false;
        }
Michael Turner's avatar
tweaks    
Michael Turner committed
336
    }
Michael Turner's avatar
Michael Turner committed
337
    else if (uv[1] > m_bounds[3])
Michael Turner's avatar
tweaks    
Michael Turner committed
338
    {
Michael Turner's avatar
Michael Turner committed
339
        if (fabs(uv[1] - m_bounds[3]) > 1E-6)
340
        {
Michael Turner's avatar
Michael Turner committed
341
342
            error << "V(" << uv[1] << ") is greater than Vmax(" << m_bounds[3]
                  << ")";
343
344
            passed = false;
        }
Michael Turner's avatar
tweaks    
Michael Turner committed
345
346
    }

347
    error << " On Surface: " << GetId();
Michael Turner's avatar
Michael Turner committed
348
    ASSERTL1(passed, "Warning: " + error.str());
Michael Turner's avatar
tweaks    
Michael Turner committed
349
}
Michael Turner's avatar
Michael Turner committed
350
}
351
}