CADSystemOCE.cpp 16.9 KB
Newer Older
mike's avatar
mike committed
1 2 3 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 32 33 34
////////////////////////////////////////////////////////////////////////////////
//
//  File: CADSystem.cpp
//
//  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.
//
//  Description: cad object methods.
//
////////////////////////////////////////////////////////////////////////////////
35
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
mike's avatar
mike committed
36

mike turner's avatar
mike turner committed
37
#include <NekMeshUtils/CADSystem/CADSurf.h>
38 39
#include <NekMeshUtils/CADSystem/OCE/CADCurveOCE.h>
#include <NekMeshUtils/CADSystem/OCE/CADSurfOCE.h>
mike turner's avatar
mike turner committed
40 41
#include <NekMeshUtils/CADSystem/OCE/CADSystemOCE.h>
#include <NekMeshUtils/CADSystem/OCE/CADVertOCE.h>
mike's avatar
mike committed
42

's avatar
committed
43
#include <boost/algorithm/string.hpp>
Michael Turner's avatar
Michael Turner committed
44
#include <boost/filesystem.hpp>
's avatar
committed
45

mike's avatar
mike committed
46 47 48 49 50 51 52
using namespace std;

namespace Nektar
{
namespace NekMeshUtils
{

53 54
std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
    "oce", CADSystemOCE::create, "Uses OCE as cad engine");
mike's avatar
mike committed
55 56 57

bool CADSystemOCE::LoadCAD()
{
58
    if (m_naca.size() == 0)
mike's avatar
mike committed
59
    {
60
        //not a naca profile behave normally
's avatar
committed
61 62 63
        //could be a geo
        string ext  = boost::filesystem::extension(m_name);

Michael Turner's avatar
Michael Turner committed
64
        if (boost::iequals(ext, ".geo"))
's avatar
committed
65 66 67 68
        {
            shape = BuildGeo(m_name);
        }
        else
69
        {
's avatar
committed
70 71 72 73 74 75 76 77 78 79 80
            // Takes step file and makes OpenCascade shape
            STEPControl_Reader reader;
            reader = STEPControl_Reader();
            reader.ReadFile(m_name.c_str());
            reader.NbRootsForTransfer();
            reader.TransferRoots();
            shape = reader.OneShape();
            if (shape.IsNull())
            {
                return false;
            }
81 82 83 84 85
        }
    }
    else
    {
        shape = BuildNACA(m_name);
mike's avatar
mike committed
86 87
    }

mike turner's avatar
mike turner committed
88
    TopExp_Explorer explr;
mike's avatar
mike committed
89 90

    // build map of verticies
mike turner's avatar
mike turner committed
91
    for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
mike's avatar
mike committed
92
    {
mike turner's avatar
mike turner committed
93
        TopoDS_Shape v = explr.Current();
94
        if (mapOfVerts.Contains(v))
mike turner's avatar
mike turner committed
95 96 97 98
        {
            continue;
        }
        int i = mapOfVerts.Add(v);
mike's avatar
mike committed
99 100 101 102 103 104
        AddVert(i, v);
    }

    // For each face of the geometry, get the local edges which bound it. If
    // they are valid (their type != 7), then add them to an edge map. This
    // filters out the dummy edges which OCC uses.
mike turner's avatar
mike turner committed
105
    for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
mike's avatar
mike committed
106
    {
mike turner's avatar
mike turner committed
107
        TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
108
        if (mapOfEdges.Contains(e))
mike's avatar
mike committed
109
        {
mike turner's avatar
mike turner committed
110
            continue;
mike's avatar
mike committed
111
        }
mike turner's avatar
mike turner committed
112 113
        BRepAdaptor_Curve curve = BRepAdaptor_Curve(TopoDS::Edge(e));
        if (curve.GetType() != 7)
mike's avatar
mike committed
114
        {
mike turner's avatar
mike turner committed
115
            int i = mapOfEdges.Add(e);
116
            AddCurve(i, e);
mike's avatar
mike committed
117 118 119
        }
    }

mike turner's avatar
mike turner committed
120
    for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
mike's avatar
mike committed
121
    {
Michael Turner's avatar
Michael Turner committed
122
        TopoDS_Shape f = explr.Current();
mike turner's avatar
mike turner committed
123 124
        ASSERTL0(!mapOfFaces.Contains(f), "duplicated faces");
        int i = mapOfFaces.Add(f);
mike's avatar
mike committed
125

mike turner's avatar
mike turner committed
126
        AddSurf(i, f);
mike's avatar
mike committed
127 128 129 130 131
    }

    // attempts to identify properties of the vertex on the degen edge
    for (int i = 1; i <= mapOfFaces.Extent(); i++)
    {
Michael Turner's avatar
Michael Turner committed
132
        TopoDS_Shape face = mapOfFaces.FindKey(i).Oriented(TopAbs_FORWARD);
mike's avatar
mike committed
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155

        TopTools_IndexedMapOfShape localEdges;
        TopExp::MapShapes(face, TopAbs_EDGE, localEdges);

        for (int j = 1; j <= localEdges.Extent(); j++)
        {
            TopoDS_Shape edge = localEdges.FindKey(j);
            if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
            {
                gp_Pnt2d p1, p2;

                BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
                                    p2);

                m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
                            TopoDS::Edge(edge), Standard_True))]
                    ->SetDegen(i, m_surfs[i], (p1.X() + p2.X()) / 2.0,
                               (p1.Y() + p2.Y()) / 2.0);
            }
        }
    }

    // This checks that all edges are bound by two surfaces, sanity check.
mike turner's avatar
mike turner committed
156
    if (!m_2d)
mike's avatar
mike committed
157
    {
mike turner's avatar
mike turner committed
158 159
        map<int, CADCurveSharedPtr>::iterator it;
        for (it = m_curves.begin(); it != m_curves.end(); it++)
160
        {
mike turner's avatar
mike turner committed
161 162
            ASSERTL0(it->second->GetAdjSurf().size() == 2,
                     "curve is not joined to 2 surfaces");
163
        }
mike's avatar
mike committed
164
    }
mike turner's avatar
mike turner committed
165

mike's avatar
mike committed
166 167 168 169 170 171 172
    return true;
}

void CADSystemOCE::AddVert(int i, TopoDS_Shape in)
{
    CADVertSharedPtr newVert = GetCADVertFactory().CreateInstance(key);

173
    boost::static_pointer_cast<CADVertOCE>(newVert)->Initialise(i, in);
mike's avatar
mike committed
174 175 176 177

    m_verts[i] = newVert;
}

mike turner's avatar
mike turner committed
178
void CADSystemOCE::AddCurve(int i, TopoDS_Shape in)
mike's avatar
mike committed
179 180
{
    CADCurveSharedPtr newCurve = GetCADCurveFactory().CreateInstance(key);
181
    boost::static_pointer_cast<CADCurveOCE>(newCurve)->Initialise(i, in);
mike's avatar
mike committed
182

mike turner's avatar
mike turner committed
183 184 185
    TopoDS_Vertex fv = TopExp::FirstVertex(TopoDS::Edge(in));
    TopoDS_Vertex lv = TopExp::LastVertex(TopoDS::Edge(in));

mike's avatar
mike committed
186
    vector<CADVertSharedPtr> vs;
mike turner's avatar
mike turner committed
187 188 189 190
    vs.push_back(m_verts[mapOfVerts.FindIndex(fv)]);
    vs.push_back(m_verts[mapOfVerts.FindIndex(lv)]);
    newCurve->SetVert(vs);

mike's avatar
mike committed
191 192 193
    m_curves[i] = newCurve;
}

mike turner's avatar
mike turner committed
194
void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
mike's avatar
mike committed
195 196
{
    CADSurfSharedPtr newSurf = GetCADSurfFactory().CreateInstance(key);
mike turner's avatar
mike turner committed
197
    boost::static_pointer_cast<CADSurfOCE>(newSurf)->Initialise(i, in);
mike's avatar
mike committed
198

199
    // do the exploration on forward oriented
Michael Turner's avatar
Michael Turner committed
200
    TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
mike turner's avatar
mike turner committed
201
    TopTools_IndexedMapOfShape mapOfWires;
Michael Turner's avatar
Michael Turner committed
202
    TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
mike turner's avatar
mike turner committed
203 204 205
    vector<EdgeLoopSharedPtr> edgeloops;
    // now we acutally analyse the loops for cad building
    for (int j = 1; j <= mapOfWires.Extent(); j++)
mike's avatar
mike committed
206
    {
mike turner's avatar
mike turner committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
        EdgeLoopSharedPtr edgeloop = EdgeLoopSharedPtr(new EdgeLoop);

        TopoDS_Shape wire = mapOfWires.FindKey(j);

        BRepTools_WireExplorer exp;

        exp.Init(TopoDS::Wire(wire));

        while (exp.More())
        {
            TopoDS_Shape edge = exp.Current();

            if (mapOfEdges.Contains(edge))
            {
                int e = mapOfEdges.FindIndex(edge);
                edgeloop->edges.push_back(m_curves[e]);
223
                edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
224 225
                                              ? CADOrientation::eForwards
                                              : CADOrientation::eBackwards);
mike turner's avatar
mike turner committed
226 227 228 229
            }
            exp.Next();
        }
        edgeloops.push_back(edgeloop);
mike's avatar
mike committed
230 231 232
    }

    int tote = 0;
mike turner's avatar
mike turner committed
233
    for (int k = 0; k < edgeloops.size(); k++)
mike's avatar
mike committed
234
    {
mike turner's avatar
mike turner committed
235
        tote += edgeloops[k]->edges.size();
mike's avatar
mike committed
236 237 238
    }

    ASSERTL0(tote != 1, "cannot handle periodic curves");
mike turner's avatar
mike turner committed
239

mike turner's avatar
mike turner committed
240 241
    CADSurf::OrientateEdges(newSurf, edgeloops);
    newSurf->SetEdges(edgeloops);
mike turner's avatar
mike turner committed
242 243

    // now the loops are orientated, tell the curves how they are
mike turner's avatar
mike turner committed
244
    for (int k = 0; k < edgeloops.size(); k++)
mike turner's avatar
mike turner committed
245
    {
mike turner's avatar
mike turner committed
246
        for (int j = 0; j < edgeloops[k]->edges.size(); j++)
mike turner's avatar
mike turner committed
247
        {
248 249
            edgeloops[k]->edges[j]->SetAdjSurf(
                make_pair(newSurf, edgeloops[k]->edgeo[j]));
mike turner's avatar
mike turner committed
250 251 252 253
        }
    }

    m_surfs[i] = newSurf;
mike's avatar
mike committed
254 255 256 257 258 259 260 261 262 263 264 265 266 267
}

Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
{
    Array<OneD, NekDouble> bound(6);
    bound[0] = numeric_limits<double>::max(); // xmin
    bound[1] = numeric_limits<double>::min(); // xmax
    bound[2] = numeric_limits<double>::max(); // ymin
    bound[3] = numeric_limits<double>::min(); // ymax
    bound[4] = numeric_limits<double>::max(); // zmin
    bound[5] = numeric_limits<double>::min(); // zmax

    for (int i = 1; i <= m_curves.size(); i++)
    {
mike turner's avatar
mike turner committed
268
        CADCurveSharedPtr c = GetCurve(i);
mike's avatar
mike committed
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283
        Array<OneD, NekDouble> ends = c->GetMinMax();

        bound[0] = min(bound[0], min(ends[0], ends[3]));
        bound[1] = max(bound[1], max(ends[0], ends[3]));

        bound[2] = min(bound[2], min(ends[1], ends[4]));
        bound[3] = max(bound[3], max(ends[1], ends[4]));

        bound[4] = min(bound[4], min(ends[2], ends[5]));
        bound[5] = max(bound[5], max(ends[2], ends[5]));
    }

    return bound;
}

284 285 286
TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
{
    ASSERTL0(naca.length() == 4, "not a 4 digit code");
287 288 289
    vector<NekDouble> data;
    ParseUtils::GenerateUnOrderedVector(m_naca.c_str(), data);
    ASSERTL0(data.size() == 5, "not a vaild domain");
290

mike turner's avatar
mike turner committed
291 292 293 294 295 296
    int n       = boost::lexical_cast<int>(naca);
    NekDouble T = (n % 100) / 100.0;
    n /= 100;
    NekDouble P = (n % 10) / 10.0;
    n /= 10;
    NekDouble M = (n % 10) / 100.0;
297

Michael Turner's avatar
working  
Michael Turner committed
298 299 300
    int np = 25;

    Array<OneD, NekDouble> xc(np);
mike turner's avatar
mike turner committed
301 302
    NekDouble dtheta = M_PI / (np - 1);
    for (int i = 0; i < np; i++)
303
    {
mike turner's avatar
mike turner committed
304
        xc[i] = (1.0 - cos(i * dtheta)) / 2.0;
305 306
    }

Michael Turner's avatar
working  
Michael Turner committed
307
    Array<OneD, NekDouble> yc(np), dyc(np);
mike turner's avatar
mike turner committed
308
    for (int i = 0; i < np; i++)
309
    {
mike turner's avatar
mike turner committed
310
        if (xc[i] < P)
311
        {
mike turner's avatar
mike turner committed
312
            yc[i]  = M / P / P * (2.0 * P * xc[i] - xc[i] * xc[i]);
313 314 315 316
            dyc[i] = 2.0 * M / P / P * (P - xc[i]);
        }
        else
        {
mike turner's avatar
mike turner committed
317 318
            yc[i] = M / (1.0 - P) / (1.0 - P) *
                    (1.0 - 2.0 * P + 2.0 * P * xc[i] - xc[i] * xc[i]);
319 320 321 322
            dyc[i] = 2.0 * M / (1.0 - P) / (1.0 - P) * (P - xc[i]);
        }
    }

Michael Turner's avatar
working  
Michael Turner committed
323
    Array<OneD, NekDouble> yt(np);
mike turner's avatar
mike turner committed
324
    for (int i = 0; i < np; i++)
325
    {
mike turner's avatar
mike turner committed
326 327 328 329
        yt[i] =
            T / 0.2 * (0.2969 * sqrt(xc[i]) - 0.1260 * xc[i] -
                       0.3516 * xc[i] * xc[i] + 0.2843 * xc[i] * xc[i] * xc[i] -
                       0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
330 331
    }

mike turner's avatar
mike turner committed
332
    Array<OneD, NekDouble> x(2 * np - 1), y(2 * np - 1);
333
    int l = 0;
mike turner's avatar
mike turner committed
334
    for (int i = np - 1; i >= 0; i--, l++)
335 336 337 338 339 340
    {
        NekDouble theta = atan(dyc[i]);

        x[l] = xc[i] - yt[i] * sin(theta);
        y[l] = yc[i] + yt[i] * cos(theta);
    }
mike turner's avatar
mike turner committed
341
    for (int i = 1; i < np; i++)
342 343 344
    {
        NekDouble theta = atan(dyc[i]);

mike turner's avatar
mike turner committed
345 346
        x[i + np - 1] = xc[i] + yt[i] * sin(theta);
        y[i + np - 1] = yc[i] - yt[i] * cos(theta);
347 348
    }

mike turner's avatar
mike turner committed
349
    TColgp_Array1OfPnt pointArray(0, 2 * np - 2);
350

mike turner's avatar
mike turner committed
351
    for (int i = 0; i < 2 * np - 1; i++)
352
    {
mike turner's avatar
mike turner committed
353
        pointArray.SetValue(i, gp_Pnt(x[i] * 1000.0, y[i] * 1000.0, 0.0));
354 355
    }

Michael Turner's avatar
working  
Michael Turner committed
356
    GeomAPI_PointsToBSpline spline(pointArray);
357 358 359 360
    Handle(Geom_BSplineCurve) curve = spline.Curve();

    BRepBuilderAPI_MakeEdge areoEdgeBuilder(curve);
    TopoDS_Edge aeroEdge = areoEdgeBuilder.Edge();
Dave Moxey's avatar
Dave Moxey committed
361
    BRepBuilderAPI_MakeEdge aeroTEBuilder(
mike turner's avatar
mike turner committed
362 363
        gp_Pnt(x[0] * 1000.0, y[0] * 1000.0, 0.0),
        gp_Pnt(x[2 * np - 2] * 1000.0, y[2 * np - 2] * 1000.0, 0.0));
364 365 366 367 368
    TopoDS_Edge TeEdge = aeroTEBuilder.Edge();

    BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
    TopoDS_Wire aeroWire = aeroWireBuilder.Wire();

369
    gp_Trsf transform;
mike turner's avatar
mike turner committed
370 371
    gp_Ax1 rotAx(gp_Pnt(500.0, 0.0, 0.0), gp_Dir(gp_Vec(0.0, 0.0, -1.0)));
    transform.SetRotation(rotAx, data[4] / 180.0 * M_PI);
372 373 374
    TopLoc_Location mv(transform);
    aeroWire.Move(mv);

mike turner's avatar
mike turner committed
375 376 377
    BRepBuilderAPI_MakeEdge domInlBuilder(
        gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0),
        gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0));
378
    TopoDS_Edge inlEdge = domInlBuilder.Edge();
379

mike turner's avatar
mike turner committed
380 381 382
    BRepBuilderAPI_MakeEdge domTopBuilder(
        gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0),
        gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0));
383
    TopoDS_Edge topEdge = domTopBuilder.Edge();
384

mike turner's avatar
mike turner committed
385 386 387
    BRepBuilderAPI_MakeEdge domOutBuilder(
        gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0),
        gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0));
388
    TopoDS_Edge outEdge = domOutBuilder.Edge();
Michael Turner's avatar
Michael Turner committed
389

mike turner's avatar
mike turner committed
390 391 392
    BRepBuilderAPI_MakeEdge domBotBuilder(
        gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0),
        gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0));
393 394 395 396 397 398 399 400
    TopoDS_Edge botEdge = domBotBuilder.Edge();

    BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
    TopoDS_Wire domWire = domWireBuilder.Wire();

    BRepBuilderAPI_MakeFace face(domWire, true);
    face.Add(aeroWire);

401 402 403
    ShapeFix_Face sf(face.Face());
    sf.FixOrientation();

404
    return sf.Face();
405 406
}

's avatar
committed
407 408 409 410 411
TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
{
    ifstream f;
    f.open(geo.c_str());

Michael Turner's avatar
Michael Turner committed
412 413 414 415 416
    map<int, string> points;
    map<int, string> lines;
    map<int, string> splines;
    map<int, string> loops;
    map<int, string> surfs;
's avatar
committed
417 418

    string fline;
Julian Marcon's avatar
Julian Marcon committed
419
    string flinetmp;
's avatar
committed
420

Michael Turner's avatar
Michael Turner committed
421
    while (!f.eof())
's avatar
committed
422
    {
Michael Turner's avatar
Michael Turner committed
423
        getline(f, fline);
's avatar
committed
424

Julian Marcon's avatar
Julian Marcon committed
425 426
        boost::erase_all(fline, "\r");

Michael Turner's avatar
Michael Turner committed
427
        if (boost::starts_with(fline, "//"))
's avatar
committed
428 429 430 431
        {
            continue;
        }

Julian Marcon's avatar
Julian Marcon committed
432
        if (!boost::contains(fline, ";"))
Julian Marcon's avatar
Julian Marcon committed
433 434 435 436 437 438 439 440
        {
            flinetmp += fline;
            continue;
        }

        fline = flinetmp + fline;
        flinetmp.clear();

's avatar
committed
441
        vector<string> tmp1, tmp2;
Michael Turner's avatar
Michael Turner committed
442
        boost::split(tmp1, fline, boost::is_any_of("="));
's avatar
committed
443

Michael Turner's avatar
Michael Turner committed
444
        boost::split(tmp2, tmp1[0], boost::is_any_of("("));
's avatar
committed
445 446

        string type = tmp2[0];
Michael Turner's avatar
Michael Turner committed
447 448
        boost::erase_all(tmp2[1], ")");
        boost::erase_all(tmp2[1], " ");
's avatar
committed
449 450
        int id = boost::lexical_cast<int>(tmp2[1]);

Michael Turner's avatar
Michael Turner committed
451 452 453 454
        boost::erase_all(tmp1[1], " ");
        boost::erase_all(tmp1[1], "{");
        boost::erase_all(tmp1[1], "}");
        boost::erase_all(tmp1[1], ";");
's avatar
committed
455 456 457

        string var = tmp1[1];

Michael Turner's avatar
Michael Turner committed
458
        if (boost::iequals(type, "Point"))
's avatar
committed
459 460 461
        {
            points[id] = var;
        }
Michael Turner's avatar
Michael Turner committed
462
        else if (boost::iequals(type, "Line"))
's avatar
committed
463 464 465
        {
            lines[id] = var;
        }
Michael Turner's avatar
Michael Turner committed
466
        else if (boost::iequals(type, "Spline"))
's avatar
committed
467 468 469
        {
            splines[id] = var;
        }
Michael Turner's avatar
Michael Turner committed
470
        else if (boost::iequals(type, "Line Loop"))
's avatar
committed
471
        {
Michael Turner's avatar
Michael Turner committed
472 473 474 475
            //line loops sometimes have negative entries for gmsh
            //orientaton purposes
            //we dont care so remove it
            boost::erase_all(var, "-");
's avatar
committed
476 477
            loops[id] = var;
        }
Michael Turner's avatar
Michael Turner committed
478
        else if (boost::iequals(type, "Plane Surface"))
's avatar
committed
479 480 481 482 483 484 485 486 487
        {
            surfs[id] = var;
        }
        else
        {
            cout << "not sure " << type << endl;
        }
    }

Michael Turner's avatar
Michael Turner committed
488
    map<int, string>::iterator it;
's avatar
committed
489

Michael Turner's avatar
Michael Turner committed
490 491 492
    // build points
    map<int, gp_Pnt> cPoints;
    for (it = points.begin(); it != points.end(); it++)
's avatar
committed
493 494 495 496
    {
        vector<NekDouble> data;
        ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);

Michael Turner's avatar
Michael Turner committed
497 498
        cPoints[it->first] =
            gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, data[2] * 1000.0);
's avatar
committed
499 500
    }

Michael Turner's avatar
Michael Turner committed
501 502 503
    // build edges
    map<int, TopoDS_Edge> cEdges;
    for (it = lines.begin(); it != lines.end(); it++)
's avatar
committed
504
    {
Michael Turner's avatar
Michael Turner committed
505
        vector<unsigned int> data;
's avatar
committed
506
        ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
Michael Turner's avatar
Michael Turner committed
507
        BRepBuilderAPI_MakeEdge em(cPoints[data[0]], cPoints[data[1]]);
's avatar
committed
508 509
        cEdges[it->first] = em.Edge();
    }
Michael Turner's avatar
Michael Turner committed
510
    for (it = splines.begin(); it != splines.end(); it++)
's avatar
committed
511
    {
Michael Turner's avatar
Michael Turner committed
512
        vector<unsigned int> data;
's avatar
committed
513 514
        ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);

Michael Turner's avatar
Michael Turner committed
515
        TColgp_Array1OfPnt pointArray(0, data.size() - 1);
's avatar
committed
516

Michael Turner's avatar
Michael Turner committed
517
        for (int i = 0; i < data.size(); i++)
's avatar
committed
518
        {
Michael Turner's avatar
Michael Turner committed
519
            pointArray.SetValue(i, cPoints[data[i]]);
's avatar
committed
520 521 522 523 524 525 526 527
        }
        GeomAPI_PointsToBSpline spline(pointArray);
        Handle(Geom_BSplineCurve) curve = spline.Curve();

        BRepBuilderAPI_MakeEdge em(curve);
        cEdges[it->first] = em.Edge();
    }

Michael Turner's avatar
Michael Turner committed
528 529 530
    // build wires
    map<int, TopoDS_Wire> cWires;
    for (it = loops.begin(); it != loops.end(); it++)
's avatar
committed
531
    {
Michael Turner's avatar
Michael Turner committed
532
        vector<unsigned int> data;
's avatar
committed
533 534
        ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
        BRepBuilderAPI_MakeWire wm;
Michael Turner's avatar
Michael Turner committed
535
        for (int i = 0; i < data.size(); i++)
's avatar
committed
536
        {
Michael Turner's avatar
Michael Turner committed
537
            wm.Add(cEdges[data[i]]);
's avatar
committed
538 539 540 541
        }
        cWires[it->first] = wm.Wire();
    }

Michael Turner's avatar
Michael Turner committed
542 543 544
    // make surface, at this point assuming its 2D (therefore only 1)
    // also going to assume that the first loop in the list is the outer domain
    ASSERTL0(surfs.size() == 1, "more than 1 surf");
's avatar
committed
545
    it = surfs.begin();
Michael Turner's avatar
Michael Turner committed
546
    vector<unsigned int> data;
's avatar
committed
547
    ParseUtils::GenerateUnOrderedVector(it->second.c_str(), data);
Michael Turner's avatar
Michael Turner committed
548
    BRepBuilderAPI_MakeFace face(cWires[data[0]], true);
Michael Turner's avatar
Michael Turner committed
549
    for (int i = 1; i < data.size(); i++)
's avatar
committed
550
    {
Michael Turner's avatar
Michael Turner committed
551
        face.Add(cWires[data[i]]);
's avatar
committed
552 553
    }

Michael Turner's avatar
Michael Turner committed
554 555
    ASSERTL0(face.Error() == BRepBuilderAPI_FaceDone, "build geo failed");

556 557
    ShapeFix_Face sf(face.Face());
    sf.FixOrientation();
558 559

    return sf.Face();
's avatar
committed
560 561
}

mike's avatar
mike committed
562 563
}
}