Branch data Line data Source code
1 : : #include "meshkit/CAMALPaver.hpp"
2 : : #include "meshkit/MeshScheme.hpp"
3 : : #include "meshkit/ModelEnt.hpp"
4 : : #include "moab/Interface.hpp"
5 : : #include "moab/ReadUtilIface.hpp"
6 : : #include "moab/Range.hpp"
7 : : #include "meshkit/RegisterMeshOp.hpp"
8 : : #include "CAMALSurfEval.hpp"
9 : : #include "CAMALSizeEval.hpp"
10 : : #include "CMLPaver.hpp"
11 : : #include <vector>
12 : :
13 : : namespace MeshKit
14 : : {
15 : :
16 : : moab::EntityType CAMALPaver::meshTps[] = {moab::MBVERTEX, moab::MBQUAD, moab::MBMAXTYPE};
17 : :
18 : 12 : CAMALPaver::CAMALPaver(MKCore *mk_core, const MEntVector &me_vec)
19 : 12 : : MeshScheme(mk_core, me_vec)
20 : : {
21 : 12 : }
22 : :
23 : 36 : CAMALPaver::~CAMALPaver()
24 : : {
25 [ - + ]: 24 : }
26 : :
27 : 13 : void CAMALPaver::setup_this()
28 : : {
29 : : // need to constrain all edges to be even
30 : 13 : constrain_even();
31 : :
32 : : // then call setup_boundary, to set up edge meshers
33 : 13 : setup_boundary();
34 : :
35 : : // then ensure that this depends on all bounding edges being meshed,
36 : : // even if the edge meshers were already set up before setup_boundary
37 : 13 : ensure_facet_dependencies(false);
38 : 13 : }
39 : :
40 : 13 : void CAMALPaver::execute_this()
41 : : {
42 : :
43 : : #ifdef HAVE_FBIGEOM
44 [ - + ]: 13 : if (mentSelection.empty())
45 : : {
46 : : // create model ents from previous op
47 : 0 : create_model_ents_from_previous_ops();
48 : : // now look at the latest SizingFunction, and set it or each model ent
49 : 0 : int latestIndexSF = 0; // maybe we would need to set it right
50 [ # # # # ]: 0 : for (MEntSelection::iterator sit = mentSelection.begin();
[ # # ]
51 : 0 : sit != mentSelection.end(); sit++) {
52 : : // make a me, for convenience
53 [ # # ]: 0 : ModelEnt *me = (*sit).first;
54 [ # # ]: 0 : me->sizing_function_index(latestIndexSF); // need to work on this one; how do we know?
55 : : }
56 : : // now, force setup of this node again, as we have added model entities to it
57 : 0 : setup_called(false);
58 : 0 : mk_core()->setup(false);
59 : : // debug
60 : 0 : mk_core()->print_graph();
61 : : // it may not be enough, we may have to execute the previous ops, that were just
62 : : // created during setup... not very clean code;
63 : 0 : mk_core()->execute_before((GraphNode *) this);
64 : : }
65 : : #endif
66 [ + - ][ + - ]: 29 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + + ]
67 : : // make a me, for convenience
68 [ + - ]: 16 : ModelEnt *me = (*sit).first;
69 : :
70 [ + - ][ + + ]: 16 : if(me->get_meshed_state()>=COMPLETE_MESH)
71 : 1 : continue;
72 : : // create a surface evaluator for this modelent, and a size evaluator
73 [ + - ]: 15 : CAMALSurfEval cse(me);
74 : :
75 [ + - ][ + - ]: 15 : SizingFunction * sz = mk_core()->sizing_function(me->sizing_function_index());
[ + - ]
76 [ + - ][ + + ]: 15 : if (!sz->variable())
77 : 13 : sz = NULL; // no function variable
78 : :
79 [ + - ][ + - ]: 30 : CAMALSizeEval mesize(me->mesh_interval_size(), sz);
80 : : // make sure the size isn't negative
81 : : // make sure the size isn't negative
82 [ + - ][ - + ]: 15 : if (mesize.get_size() == -1) mesize.set_size(1.0);
[ # # ]
83 : :
84 : : // assemble bounding mesh
85 [ + - ]: 30 : std::vector<moab::EntityHandle> bdy;
86 [ + - ][ + - ]: 30 : std::vector<int> group_sizes, bdy_ids;
87 [ + - ]: 15 : me->boundary(0, bdy, NULL, &group_sizes);
88 : :
89 : : // convert the handles to integers for input to Paver
90 [ + - ]: 30 : moab::Range bdy_vrange;
91 [ + - ]: 30 : std::vector<double> coords;
92 [ + - ]: 15 : me->get_indexed_connect_coords(bdy, NULL, NULL, bdy_ids, coords, &bdy_vrange);
93 : :
94 : : // now construct the CAMAL mesher, and pass it initial conditions
95 [ + - ]: 30 : CMLPaver paver(&cse, &mesize);
96 [ + - ][ + - ]: 15 : bool success = paver.set_boundary_mesh(bdy_vrange.size(), &coords[0], group_sizes.size(), &group_sizes[0], &bdy_ids[0]);
[ + - ][ + - ]
[ + - ]
97 [ - + ]: 15 : if (!success)
98 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Trouble setting boundary mesh.");
99 : :
100 : : // ok, now generate the mesh
101 : : int num_pts, num_quads;
102 [ + - ]: 15 : success = paver.generate_mesh(num_pts, num_quads);
103 [ - + ]: 15 : if (!success)
104 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Trouble generating quad mesh.");
105 : :
106 [ + - ]: 15 : moab::Range &new_ents = (*sit).second;
107 : : moab::ErrorCode rval;
108 : :
109 : : // resize the coords array, then get the coords of the new points
110 [ + - ][ - + ]: 15 : assert(num_pts >= (int)bdy_vrange.size());
111 [ + - ][ + - ]: 15 : if (num_pts > (int)bdy_vrange.size()) {
112 [ + - ][ + - ]: 15 : coords.resize(3*(num_pts-bdy_vrange.size()));
113 [ + - ][ + - ]: 15 : int pts_returned = paver.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
[ + - ]
114 [ + - ][ - + ]: 15 : if (pts_returned != num_pts-(int)bdy_vrange.size())
115 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Number of new points returned from Paver doesn't agree with previous value output.");
116 : :
117 : : // create the new vertices' entities on the face
118 [ + - ][ + - ]: 15 : rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
[ + - ][ + - ]
119 [ - + ][ # # ]: 15 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
120 : : }
121 : :
122 : : // for quads, pre-allocate connectivity
123 : : moab::ReadUtilIface *iface;
124 [ + - ][ + - ]: 15 : rval = mk_core()-> moab_instance() -> query_interface(iface);
[ + - ]
125 [ - + ][ # # ]: 15 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
126 : :
127 : : //create the quads, get a direct ptr to connectivity
128 : : moab::EntityHandle starth, *connect;
129 [ + - ]: 15 : rval = iface->get_element_connect(num_quads, 4, moab::MBQUAD, 1, starth, connect);
130 [ - + ][ # # ]: 15 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
131 : :
132 : : // read connectivity directly into that array, as int's
133 : 15 : int *connecti = (int*) connect;
134 [ + - ]: 15 : int quads_returned = paver.get_quads_buf(4*num_quads, connecti);
135 [ - + ]: 15 : if (quads_returned != num_quads)
136 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Number of new quads returned from Paver doesn't agree with previous value output.");
137 : :
138 : : // put vertex handles into an indexible array
139 [ + - ]: 30 : std::vector<moab::EntityHandle> bdy_vvec;
140 [ + - ][ + - ]: 15 : std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
141 [ + - ][ + - ]: 15 : std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
142 : :
143 : : // now convert vertex indices into handles in-place, working from the back
144 [ + + ]: 6599 : for (int i = 4*num_quads-1; i >= 0; i--) {
145 [ + - ][ - + ]: 6584 : assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
146 [ + - ]: 6584 : connect[i] = bdy_vvec[connecti[i]];
147 : : }
148 : :
149 : : // put new quads into new entity range, then commit the mesh
150 [ + - ]: 15 : new_ents.insert(starth, starth+num_quads-1);
151 [ + - ]: 15 : me->commit_mesh(new_ents, COMPLETE_MESH);
152 : 15 : }
153 : 13 : }
154 : :
155 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
