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
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#include "meshkit/TGTetMesher.hpp"
#include "meshkit/MeshScheme.hpp"
#include "meshkit/ModelEnt.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
#include "meshkit/RegisterMeshOp.hpp"

//namespace tetgen {
  #include "tetgen.h"
//}



#include <vector>

namespace MeshKit
{

moab::EntityType TGTetMesher::meshTps[] = {moab::MBVERTEX, moab::MBTET, moab::MBMAXTYPE};

TGTetMesher::TGTetMesher(MKCore *mk_core, const MEntVector &me_vec)
        : MeshScheme(mk_core, me_vec)
{
  tet_switch = (char*)"pq1.414Y";
}

TGTetMesher::~TGTetMesher()
{
}

MeshOp *TGTetMesher::get_tri_mesher()
{
  std::vector<MeshOpProxy *> proxies;<--- Unused variable: proxies
  //  mk_core()->meshop_by_mesh_type(moab::MBTRI, proxies);
  // if (proxies.empty()) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing triangles.");
  //return mk_core()->construct_meshop(*proxies.begin());
  return mk_core()->construct_meshop("CAMALTriAdvance");
}

void TGTetMesher::setup_this()
{
  MeshOp *tri_mesher = NULL;
  MEntVector surfs;
  
  for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
    ModelEnt *me = (*sit).first;
    if (me->dimension() != 3) throw Error(MK_BAD_INPUT, "Tet mesher assigned to an entity with dimension != 3.");
    
      // get the bounding faces
    surfs.clear();
    me->get_adjacencies(2, surfs);
    
      // check the mesh scheme on each one; if there is one, verify it can generate tris; if there
      // isn't one, make one
    bool inserted = false;
    
    for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
      if ((*vit)->is_meshops_list_empty()) {
          // get a tri mesher if we haven't already
        if (!tri_mesher) tri_mesher = get_tri_mesher();
        
          // add this surface to it, and if first for the volume, make sure it's added upstream
        tri_mesher->add_modelent(*vit);
        if (!inserted) {
          mk_core()->insert_node(tri_mesher, this);
          inserted = true;
        }
      } // if no meshops
    } // over surfs
  } // over vols
}

void TGTetMesher::execute_this()
{
  tetgenio in, out;
  tetgenio::facet *f;
  tetgenio::polygon *p;
  in.firstnumber = 0;

  for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
      // make a me, for convenience
    ModelEnt *me = (*sit).first;

    // assemble bounding mesh
    std::vector<moab::EntityHandle> bdy;
    std::vector<int> group_sizes, senses, bdy_ids;
    me->boundary(2, bdy, &senses, &group_sizes);
    
    // get connectivity

      // convert the handles to integers for input to CMLTetMesher
    moab::Range bdy_vrange;
    std::vector<double> coords;
    me->get_indexed_connect_coords(bdy, &senses, NULL, bdy_ids, coords, &bdy_vrange, 1);

    // add points
    in.numberofpoints = bdy_vrange.size();
    in.pointlist = new REAL[3 * bdy_vrange.size()];
    for (int i = 0; i < (int)bdy_vrange.size(); i++)
    {
        in.pointlist[3 * i + 0] = coords[3*i];
        in.pointlist[3 * i + 1] = coords[3*i+1];
        in.pointlist[3 * i + 2] = coords[3*i+2];
      }

     // add the tris
     unsigned int numtris = bdy.size();

     in.numberoffacets = numtris;
     in.facetlist = new tetgenio::facet[numtris];
     in.facetmarkerlist = NULL;

    for (unsigned int i = 0; i < numtris; i++){
        f = &in.facetlist[i];
        f->numberofpolygons = 1;
        f->polygonlist = new tetgenio::polygon[1];
        f->numberofholes = 0;
        f->holelist = NULL;
        p = &f->polygonlist[0];
        p->numberofvertices = 3;
        p->vertexlist = new int[3];
        p->vertexlist[0] = bdy_ids[3*i];
        p->vertexlist[1] = bdy_ids[3*i + 1];
        p->vertexlist[2] = bdy_ids[3*i + 2];
      }

    // CALL TETGEN
    tetrahedralize(tet_switch, &in, &out);

   // TODO: Bring back mesh to MeshKit database
    }
}

} // namespace MeshKit