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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155 | #include "meshkit/CAMALPaver.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"
#include "CAMALSurfEval.hpp"
#include "CAMALSizeEval.hpp"
#include "CMLPaver.hpp"
#include <vector>
namespace MeshKit
{
moab::EntityType CAMALPaver::meshTps[] = {moab::MBVERTEX, moab::MBQUAD, moab::MBMAXTYPE};
CAMALPaver::CAMALPaver(MKCore *mk_core, const MEntVector &me_vec)
: MeshScheme(mk_core, me_vec)
{
}
CAMALPaver::~CAMALPaver()
{
}
void CAMALPaver::setup_this()
{
// need to constrain all edges to be even
constrain_even();
// then call setup_boundary, to set up edge meshers
setup_boundary();
// then ensure that this depends on all bounding edges being meshed,
// even if the edge meshers were already set up before setup_boundary
ensure_facet_dependencies(false);
}
void CAMALPaver::execute_this()
{
#ifdef HAVE_FBIGEOM
if (mentSelection.empty())
{
// create model ents from previous op
create_model_ents_from_previous_ops();
// now look at the latest SizingFunction, and set it or each model ent
int latestIndexSF = 0; // maybe we would need to set it right
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;
me->sizing_function_index(latestIndexSF); // need to work on this one; how do we know?
}
// now, force setup of this node again, as we have added model entities to it
setup_called(false);
mk_core()->setup(false);
// debug
mk_core()->print_graph();
// it may not be enough, we may have to execute the previous ops, that were just
// created during setup... not very clean code;
mk_core()->execute_before((GraphNode *) this);
}
#endif
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;
if(me->get_meshed_state()>=COMPLETE_MESH)
continue;
// create a surface evaluator for this modelent, and a size evaluator
CAMALSurfEval cse(me);
SizingFunction * sz = mk_core()->sizing_function(me->sizing_function_index());
if (!sz->variable())
sz = NULL; // no function variable
CAMALSizeEval mesize(me->mesh_interval_size(), sz);
// make sure the size isn't negative
// make sure the size isn't negative
if (mesize.get_size() == -1) mesize.set_size(1.0);
// assemble bounding mesh
std::vector<moab::EntityHandle> bdy;
std::vector<int> group_sizes, bdy_ids;
me->boundary(0, bdy, NULL, &group_sizes);
// convert the handles to integers for input to Paver
moab::Range bdy_vrange;
std::vector<double> coords;
me->get_indexed_connect_coords(bdy, NULL, NULL, bdy_ids, coords, &bdy_vrange);
// now construct the CAMAL mesher, and pass it initial conditions
CMLPaver paver(&cse, &mesize);
bool success = paver.set_boundary_mesh(bdy_vrange.size(), &coords[0], group_sizes.size(), &group_sizes[0], &bdy_ids[0]);
if (!success)
ECERRCHK(MK_FAILURE, "Trouble setting boundary mesh.");
// ok, now generate the mesh
int num_pts, num_quads;
success = paver.generate_mesh(num_pts, num_quads);
if (!success)
ECERRCHK(MK_FAILURE, "Trouble generating quad mesh.");
moab::Range &new_ents = (*sit).second;
moab::ErrorCode rval;
// resize the coords array, then get the coords of the new points
assert(num_pts >= (int)bdy_vrange.size());
if (num_pts > (int)bdy_vrange.size()) {
coords.resize(3*(num_pts-bdy_vrange.size()));
int pts_returned = paver.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
if (pts_returned != num_pts-(int)bdy_vrange.size())
ECERRCHK(MK_FAILURE, "Number of new points returned from Paver doesn't agree with previous value output.");
// create the new vertices' entities on the face
rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
MBERRCHK(rval, mk_core()->moab_instance());
}
// for quads, pre-allocate connectivity
moab::ReadUtilIface *iface;
rval = mk_core()-> moab_instance() -> query_interface(iface);
MBERRCHK(rval, mk_core()->moab_instance());
//create the quads, get a direct ptr to connectivity
moab::EntityHandle starth, *connect;
rval = iface->get_element_connect(num_quads, 4, moab::MBQUAD, 1, starth, connect);
MBERRCHK(rval, mk_core()->moab_instance());
// read connectivity directly into that array, as int's
int *connecti = (int*) connect;
int quads_returned = paver.get_quads_buf(4*num_quads, connecti);
if (quads_returned != num_quads)
ECERRCHK(MK_FAILURE, "Number of new quads returned from Paver doesn't agree with previous value output.");
// put vertex handles into an indexible array
std::vector<moab::EntityHandle> bdy_vvec;
std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
// now convert vertex indices into handles in-place, working from the back
for (int i = 4*num_quads-1; i >= 0; i--) {
assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
connect[i] = bdy_vvec[connecti[i]];
}
// put new quads into new entity range, then commit the mesh
new_ents.insert(starth, starth+num_quads-1);
me->commit_mesh(new_ents, COMPLETE_MESH);
}
}
} // namespace MeshKit
|