Branch data Line data Source code
1 : : #include "meshkit/iGeom.hpp"
2 : : #include "meshkit/ModelEnt.hpp"
3 : : #include "meshkit/SizingFunction.hpp"
4 : : #include "meshkit/Error.hpp"
5 : : #include "meshkit/MKCore.hpp"
6 : : #include "meshkit/Types.hpp"
7 : : #include "meshkit/VecUtil.hpp"
8 : : #include "moab/GeomTopoTool.hpp"
9 : : #include "moab/CN.hpp"
10 : : #include "RefEntity.hpp"
11 : :
12 : : namespace MeshKit
13 : : {
14 : :
15 : 1188 : ModelEnt::ModelEnt(MKCore *mk,
16 : : iGeom::EntityHandle geom_ent,
17 : : int geom_index,
18 : : moab::EntityHandle mesh_ent,
19 : : int mesh_index,
20 : : int irel_index,
21 : : int sizing_index)
22 : : : mkCore(mk), igeomIndex(geom_index), meshIndex(mesh_index), irelIndex(irel_index),
23 : : iGeomEnt(geom_ent), iGeomSet(NULL),
24 : : moabEntSet(mesh_ent), sizingFunctionIndex(sizing_index),
25 : : meshIntervals(-1), intervalFirmness(DEFAULT), constrainEven(false),
26 : 1188 : meshedState(NO_MESH), iaVariable(NULL)
27 : : {
28 : : // mark the geometry entity with this modelEnt
29 [ + - ][ + + ]: 1188 : if (-1 != igeomIndex && geom_ent) {
30 : 1141 : ModelEnt *dum_this = this;
31 [ + - ][ + - ]: 1141 : iGeom::Error err = igeom_instance()->setData(geom_ent, mkCore->igeom_model_tag(igeomIndex), &dum_this);
[ + - ]
32 [ + - ]: 1141 : IBERRCHK(err, "Failed to set iGeom ModelEnt tag.");
33 : : }
34 : :
35 [ + + ][ + - ]: 1188 : if (!mesh_ent && geom_ent && -1 != meshIndex && -1 != irelIndex) {
[ + - ][ + + ]
36 : : // get a corresponding mesh set handle, if any; create one if missing and have a meshIndex
37 : : iBase_EntitySetHandle h;
38 [ + - ][ + - ]: 1096 : iRel::Error err = mkCore->irel_pair(irelIndex)->getEntSetRelation(geom_ent, false, h);
39 [ + + ]: 1096 : if (iBase_SUCCESS == err)
40 : 144 : moabEntSet = reinterpret_cast<moab::EntityHandle>(h);
41 : 1096 : else moabEntSet = 0;
42 : : }
43 : :
44 : : // if still no mesh entity and non-default mesh index, create one
45 [ + - ][ + + ]: 1188 : if (-1 != meshIndex && !moabEntSet) {
46 : : // should have a geometry entity here
47 [ - + ][ # # ]: 997 : assert(iGeomEnt || iGeomSet);
48 [ + - ]: 997 : create_mesh_set();
49 : : }
50 : :
51 : : // if there's a mesh entity, set it to point here too
52 [ + - ][ + - ]: 1188 : if (moabEntSet && -1 != meshIndex) {
53 : 1188 : ModelEnt *dum_this = this;
54 [ + - ]: 1188 : moab::ErrorCode err = mkCore->moab_instance(meshIndex)->tag_set_data(mkCore->moab_model_tag(meshIndex),
55 [ + - ][ + - ]: 1188 : &moabEntSet, 1, &dum_this);
56 [ + - ]: 1188 : IBERRCHK(err, "Failed to set moab ModelEnt tag.");
57 : : }
58 : 1188 : }
59 : :
60 : 6 : ModelEnt::ModelEnt(MKCore *mk,
61 : : iGeom::EntitySetHandle geom_ent,
62 : : int geom_index,
63 : : moab::EntityHandle mesh_ent,
64 : : int mesh_index,
65 : : int irel_index,
66 : : int sizing_index,
67 : : IAVariable *ia_var)
68 : : : mkCore(mk), igeomIndex(geom_index), meshIndex(mesh_index), irelIndex(irel_index),
69 : : iGeomEnt(NULL), iGeomSet(geom_ent),
70 : : moabEntSet(mesh_ent), sizingFunctionIndex(sizing_index),
71 : : meshIntervals(-1), intervalFirmness(DEFAULT), constrainEven(false),
72 : 6 : meshedState(NO_MESH), iaVariable(ia_var)
73 : : {
74 : : // mark the geometry entity with this modelEnt
75 [ + - ][ - + ]: 6 : if (-1 != igeomIndex && geom_ent) {
76 : 0 : ModelEnt *dum_this = this;
77 [ # # ]: 0 : iGeom::Error err = igeom_instance()->setEntSetData(geom_ent,
78 : : mkCore->igeom_model_tag(igeomIndex),
79 [ # # ][ # # ]: 0 : &dum_this);
80 [ # # ]: 0 : IBERRCHK(err, "Failed to set iGeom ModelSet tag.");
81 : : }
82 : :
83 [ + + ][ - + ]: 6 : if (!mesh_ent && geom_ent && -1 != meshIndex && -1 != irelIndex) {
[ # # ][ # # ]
84 : : // get a corresponding mesh set handle, if any; create one if missing and have a meshIndex
85 : : iBase_EntitySetHandle h;
86 [ # # ][ # # ]: 0 : iRel::Error err = mkCore->group_set_pair(irelIndex)->getSetSetRelation(geom_ent, false, h);
87 : : //IBERRCHK(err, "Failed to get mesh set handle for geometry entity set.");
88 [ # # ]: 0 : if(iBase_SUCCESS == err) moabEntSet = reinterpret_cast<moab::EntityHandle>(h);
89 : 0 : else moabEntSet=0;
90 : : }
91 : :
92 : : // if still no mesh entity and non-default mesh index, create one
93 [ + + ][ - + ]: 6 : if (-1 != meshIndex && !moabEntSet) {
94 : : // should have a geometry entity here
95 [ # # ][ # # ]: 0 : assert(iGeomEnt || iGeomSet);
96 [ # # ]: 0 : create_mesh_set();
97 : : }
98 : :
99 : : // if there's a mesh entity, set it to point here too
100 [ + + ][ + + ]: 6 : if (moabEntSet && -1 != meshIndex) {
101 : 2 : ModelEnt *dum_this = this;
102 [ + - ]: 2 : moab::ErrorCode err = mkCore->moab_instance(meshIndex)->tag_set_data(mkCore->moab_model_tag(meshIndex),
103 [ + - ][ + - ]: 2 : &moabEntSet, 1, &dum_this);
104 [ + - ]: 2 : IBERRCHK(err, "Failed to set moab ModelEnt tag.");
105 : : }
106 : 6 : }
107 : :
108 : 1925 : ModelEnt::~ModelEnt()
109 : : {
110 : : // if we delete the Model Entity, remove the tag on the gent that
111 : : // was pointing to it (or set it to NULL)
112 : : /*if (-1 != igeomIndex && iGeomEnt) {
113 : : iGeom::Error err = igeom_instance()->rmvTag(iGeomEnt,
114 : : mkCore->igeom_model_tag(igeomIndex));
115 : : IBERRCHK(err, "Failed to set iGeom ModelEnt tag.");
116 : : }
117 : :
118 : : if (-1 != igeomIndex && iGeomSet) {
119 : : iGeom::Error err = igeom_instance()->rmvEntSetTag(iGeomSet,
120 : : mkCore->igeom_model_tag(igeomIndex));
121 : : IBERRCHK(err, "Failed to set iGeom ModelSet tag.");
122 : : }*/
123 : : // if we delete the Model Entity, remove the tag on the moabEntSet
124 : : // that was pointing to it (set the value to NULL)
125 [ + + ][ + + ]: 643 : if (moabEntSet && -1 != meshIndex) {
126 : 639 : moab::ErrorCode err = mkCore->moab_instance(meshIndex)->tag_delete_data(
127 : 639 : mkCore->moab_model_tag(meshIndex), &moabEntSet, 1);
128 : 639 : IBERRCHK(err, "Failed to set moab ModelEnt tag.");
129 : : }
130 [ - + ]: 1282 : }
131 : :
132 : 597 : void ModelEnt::sizing_function_index(int ind, bool children_too)
133 : : {
134 : 597 : sizingFunctionIndex = ind;
135 : :
136 : : // set on children too, if requested
137 [ + + ][ + + ]: 597 : if (children_too && dimension() > 1) {
[ + + ]
138 [ + - ]: 102 : MEntVector childrn;
139 [ + - ][ + - ]: 102 : get_adjacencies(dimension()-1, childrn);
140 [ + - ][ + - ]: 622 : for (MEntVector::iterator vit = childrn.begin(); vit != childrn.end(); vit++)
[ + + ]
141 [ + - ][ + - ]: 622 : (*vit)->sizing_function_index(ind, children_too);
142 : : }
143 : 597 : }
144 : :
145 : 997 : void ModelEnt::create_mesh_set(int ordered_flag)
146 : : {
147 [ - + ]: 997 : if (moabEntSet)
148 [ # # ]: 0 : throw Error(MK_FAILURE, "Tried to create meshset for an entity that already had one.");
149 : :
150 [ - + ][ # # ]: 997 : else if (!iGeomEnt && !iGeomSet)
151 [ # # ]: 0 : throw Error(MK_FAILURE, "Tried to create ModelEnt missing a geometry entity or set.");
152 : :
153 : 997 : moab::EntitySetProperty ordered = moab::MESHSET_SET;
154 : :
155 : : // need type for later
156 : 997 : iBase_EntityType this_tp = iBase_ALL_TYPES;
157 : : iGeom::Error err;
158 [ + - ][ - + ]: 997 : if (-1 > ordered_flag || 1 < ordered_flag) throw Error(MK_FAILURE, "Invalid value for ordered flag.");
[ # # ]
159 [ - + ]: 997 : else if (iGeomSet) ordered = moab::MESHSET_ORDERED;
160 [ + - ][ + - ]: 997 : else if (-1 == ordered_flag && iGeomEnt) {
161 [ + - ][ + - ]: 997 : err = igeom_instance()->getEntType(iGeomEnt, this_tp);
162 [ + - ]: 997 : IBERRCHK(err, "Trouble getting entity type.");
163 [ + + ]: 997 : if (iBase_EDGE == this_tp)
164 : 449 : ordered = moab::MESHSET_ORDERED;
165 : : }
166 : :
167 : : // create the set
168 [ + - ][ + - ]: 997 : moab::ErrorCode rval = mkCore->moab_instance(meshIndex)->create_meshset(ordered, moabEntSet);
169 [ - + ][ # # ]: 997 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
170 : :
171 : : // set dimension tag and global id tag
172 [ + - ]: 997 : if (iBase_ALL_TYPES != this_tp) {
173 [ + - ][ + - ]: 997 : rval = mkCore->moab_instance(meshIndex)->tag_set_data(mkCore->moab_geom_dim_tag(), &moabEntSet, 1, &this_tp);
[ + - ]
174 [ - + ][ # # ]: 997 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
175 : : // set the id tag; get it with id(); if it returns -1, set a value 0
176 [ + - ]: 997 : int id_from_geometry = id();
177 [ - + ]: 997 : if (id_from_geometry==-1)
178 : 0 : id_from_geometry = 0;
179 [ + - ][ + - ]: 997 : rval = mkCore->moab_instance(meshIndex)->tag_set_data(mkCore->moab_global_id_tag(), &moabEntSet, 1, &id_from_geometry);
[ + - ]
180 [ - + ][ # # ]: 997 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
181 : : }
182 : :
183 : : // tag it with this ModelEnt
184 : 997 : ModelEnt *this_ptr = this;
185 [ + - ][ + - ]: 997 : rval = mkCore->moab_instance()->tag_set_data(mkCore->moab_model_tag(), &moabEntSet, 1, &this_ptr);
[ + - ]
186 [ - + ][ # # ]: 997 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
187 : :
188 : : // tag it with geometry dimension
189 : 997 : int geom_dim = this_tp;
190 [ + - ][ + - ]: 997 : rval = mkCore->moab_instance()->tag_set_data(mkCore->moab_geom_dim_tag(), &moabEntSet, 1, &geom_dim);
[ + - ]
191 [ - + ][ # # ]: 997 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
192 : :
193 : : // relate the mesh to the geom, only if iRelFlag is true
194 [ + + ]: 997 : if (-1 != irelIndex)
195 : : {
196 : : iRel::Error merr;
197 [ + - ]: 952 : if (iGeomEnt) {
198 [ + - ][ + - ]: 952 : merr = mkCore->irel_pair(irelIndex)->setEntSetRelation(iGeomEnt, IBSH(moabEntSet));
199 : : }
200 : : else {
201 [ # # ][ # # ]: 0 : merr = mkCore->group_set_pair()->setSetSetRelation(iGeomSet, IBSH(moabEntSet));
202 : : }
203 [ + - ]: 952 : IBERRCHK(merr, "Failed to set iRel relation for a mesh set.");
204 : : }
205 : :
206 [ + - ][ + - ]: 997 : if (iGeomEnt) init_parents_children();
207 [ # # ][ # # ]: 0 : else if (iGeomSet) init_group_contents();
208 : :
209 : : // now do senses wrt lower-dimensional entities, which should be initialized by now
210 [ + - ]: 997 : if (iGeomEnt)
211 [ + - ]: 997 : set_downward_senses();
212 : 997 : }
213 : :
214 : 997 : void ModelEnt::init_parents_children()
215 : : {
216 : : // get parents and children, and link to corresponding sets; don't do this for any missing mesh sets,
217 : : // will get done when those sets get created
218 [ - + ]: 997 : assert(iGeomEnt);
219 [ + - ]: 997 : std::vector<iGeom::EntityHandle> geom_adjs;
220 : 997 : std::vector<iGeom::EntityHandle>::iterator vit;
221 : : moab::EntityHandle mseth;
222 : : //moab::GeomTopoTool gt(mkCore->moab_instance(meshIndex), false);
223 : : iBase_EntityType this_tp;
224 : : moab::ErrorCode rval;
225 [ + - ][ + - ]: 997 : iGeom::Error err = igeom_instance()->getEntType(iGeomEnt, this_tp);
226 [ + - ]: 997 : IBERRCHK(err, "Trouble getting entity type.");
227 [ + + ]: 997 : if (this_tp < iBase_REGION) {
228 [ + - ][ + - ]: 956 : err = igeom_instance()->getEntAdj(iGeomEnt, (iBase_EntityType)(this_tp+1), geom_adjs);
229 [ + - ]: 956 : IBERRCHK(err, "Trouble finding parent entities.");
230 [ + - ][ + - ]: 2817 : for (vit = geom_adjs.begin(); vit != geom_adjs.end(); vit++) {
[ + + ]
231 [ + - ][ + + ]: 1861 : try {mseth = mesh_handle(*vit);
232 : : }
233 : 3590 : catch (Error err) {
234 : : // just continue, will get this when parent gets created later
235 : 1795 : continue;
236 : : }
237 [ - + ]: 66 : if (mseth) {
238 [ # # ][ # # ]: 0 : rval = mkCore->moab_instance(meshIndex)->add_parent_child(mseth, moabEntSet);
239 [ # # ][ # # ]: 0 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
240 : : }
241 : : }
242 : : }
243 : :
244 [ + + ]: 997 : if (this_tp > iBase_VERTEX) {
245 : 667 : geom_adjs.clear();
246 [ + - ][ + - ]: 667 : err = igeom_instance()->getEntAdj(iGeomEnt, (iBase_EntityType)(this_tp-1), geom_adjs);
247 [ + - ]: 667 : IBERRCHK(err, "Trouble finding child entities.");
248 [ + - ][ + - ]: 2528 : for (vit = geom_adjs.begin(); vit != geom_adjs.end(); vit++) {
[ + + ]
249 : : // should definitely have a mesh handle here, error if not
250 [ + - ][ + - ]: 1861 : mseth = mesh_handle(*vit);
251 [ + - ]: 1861 : if (mseth) {
252 [ + - ][ + - ]: 1861 : rval = mkCore->moab_instance(meshIndex)->add_parent_child(moabEntSet, mseth);
253 [ - + ][ # # ]: 1861 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
254 : : }
255 : : }
256 : 997 : }
257 [ - + ]: 2792 : }
258 : :
259 : 0 : void ModelEnt::init_group_contents()
260 : : {
261 : : // should only be calling this function if we have a geometry set and a mesh index
262 [ # # ][ # # ]: 0 : assert(-1 != meshIndex && iGeomSet);
263 : :
264 : : // get the geometry entities in this group and add the corresponding mesh entity sets
265 [ # # ]: 0 : std::vector<iGeom::EntityHandle> geom_adjs;
266 : 0 : std::vector<iGeom::EntityHandle>::iterator vit;
267 : : moab::EntityHandle mseth;
268 : : moab::ErrorCode rval;
269 [ # # ][ # # ]: 0 : iGeom::Error err = igeom_instance()->getEntities(iGeomSet, iBase_ALL_TYPES, geom_adjs);
270 [ # # ]: 0 : IBERRCHK(err, "Trouble finding parent entities.");
271 [ # # ][ # # ]: 0 : for (vit = geom_adjs.begin(); vit != geom_adjs.end(); vit++) {
[ # # # # ]
272 : 0 : mseth = 0;
273 [ # # ][ # # ]: 0 : try {mseth = mesh_handle(*vit);
274 : : }
275 : 0 : catch (Error) {
276 : : // just continue, will get this when parent gets created later
277 : 0 : continue;
278 : : }
279 [ # # ]: 0 : if (mseth) {
280 [ # # ][ # # ]: 0 : rval = mkCore->moab_instance(meshIndex)->add_entities(moabEntSet, &mseth, 1);
281 [ # # ][ # # ]: 0 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
282 : : }
283 : : }
284 : :
285 : : // same for set members
286 [ # # ]: 0 : std::vector<iGeom::EntitySetHandle> geom_sadjs;
287 : 0 : std::vector<iGeom::EntitySetHandle>::iterator svit;
288 [ # # ][ # # ]: 0 : err = igeom_instance()->getEntSets(iGeomSet, -1, geom_sadjs);
289 [ # # ]: 0 : IBERRCHK(err, "Trouble finding parent entities.");
290 [ # # ][ # # ]: 0 : for (svit = geom_sadjs.begin(); svit != geom_sadjs.end(); svit++) {
[ # # ]
291 : 0 : mseth = 0;
292 [ # # ][ # # ]: 0 : try {mseth = mesh_handle(*svit);
293 : : }
294 : 0 : catch (Error) {
295 : : // just continue, will get this when parent gets created later
296 : 0 : continue;
297 : : }
298 [ # # ]: 0 : if (mseth) {
299 [ # # ][ # # ]: 0 : rval = mkCore->moab_instance(meshIndex)->add_entities(moabEntSet, &mseth, 1);
300 [ # # ][ # # ]: 0 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
301 : : }
302 : 0 : }
303 [ # # ]: 0 : }
304 : :
305 : 997 : void ModelEnt::set_downward_senses()
306 : : {
307 : : // should have both mesh and geometry
308 [ + - ][ + - ]: 997 : assert(iGeomEnt && moabEntSet && igeomIndex != -1 && meshIndex != -1);
[ + - ][ - + ]
309 : :
310 : : // skip if vertex or edge
311 [ + - ][ + + ]: 1215 : if (dimension() <= 1) return;
312 : :
313 [ + - ]: 218 : int dim = dimension();
314 : : iGeom::Error err;
315 [ + - ]: 218 : MEntVector adjs;
316 [ + - ]: 218 : get_adjacencies(dim-1, adjs);
317 : :
318 [ + - ]: 436 : std::vector<int> senses;
319 : : int dum_sense;
320 [ + - ]: 436 : std::vector<moab::EntityHandle> ments;
321 [ + - ][ + - ]: 436 : moab::GeomTopoTool gt(mkCore->moab_instance(meshIndex));
322 [ + - ][ + - ]: 1199 : for (MEntVector::iterator vit = adjs.begin(); vit != adjs.end(); vit++) {
[ + + ]
323 [ + - ][ + - ]: 981 : err = igeom_instance()->getSense((*vit)->geom_handle(), geom_handle(), dum_sense);
[ + - ][ + - ]
[ + - ]
324 [ + - ]: 981 : IBERRCHK(err, "Problem getting senses.");
325 [ + - ][ + - ]: 981 : moab::ErrorCode rval = gt.set_sense((*vit)->mesh_handle(), mesh_handle(), dum_sense);
[ + - ][ + - ]
326 [ - + ][ # # ]: 981 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
327 : 218 : }
328 : : }
329 : :
330 : : /** \brief Commit mesh to a model entity
331 : : *
332 : : * Takes the input mesh entities, adds them to the entity set for this model entity,
333 : : * and (if both-type relation on the mesh side) sets the relations to the corresponding
334 : : * geometry entity.
335 : : * \param mesh_ents Mesh entities being assigned to this model entity
336 : : * \param mstate The meshed state after this mesh is added
337 : : */
338 : 631 : void ModelEnt::commit_mesh(moab::Range &mesh_ents, MeshedState mstate)
339 : : {
340 [ + - ]: 631 : std::vector<moab::EntityHandle> ent_vec;
341 [ + - ][ + - ]: 631 : std::copy(mesh_ents.begin(), mesh_ents.end(), std::back_inserter(ent_vec));
[ + - ][ + - ]
342 [ + - ][ + - ]: 631 : commit_mesh(&ent_vec[0], ent_vec.size(), mstate);
343 : 631 : }
344 : :
345 : : /** \brief Commit mesh to a model entity
346 : : *
347 : : * Takes the input mesh entities, and adds them to the entity set for this model entity.
348 : : * \param mesh_ents Mesh entities being assigned to this model entity
349 : : * \param mstate The meshed state after this mesh is added
350 : : */
351 : 647 : void ModelEnt::commit_mesh(moab::EntityHandle *mesh_ents,
352 : : int num_ents,
353 : : MeshedState mstate)
354 : : {
355 [ - + ]: 647 : assert(-1 != meshIndex);
356 : :
357 : : // add them to the set
358 : 647 : moab::ErrorCode rval = mkCore->moab_instance(meshIndex)->add_entities(moabEntSet, mesh_ents, num_ents);
359 [ - + ][ # # ]: 647 : MBERRCHK(rval, mkCore->moab_instance(meshIndex));
[ # # ][ # # ]
[ # # ]
360 : :
361 : 647 : meshedState = mstate;
362 : 647 : }
363 : :
364 : 4156 : void ModelEnt::get_adjacencies(int dim, MEntVector &adjs) const
365 : : {
366 [ + - ]: 4156 : std::vector<iGeom::EntityHandle> gents;
367 : : // this call assumes that we have geometry. Is that always true?
368 : : // it would be possible to have ModelEnt without geometry (geometry handle null),
369 : : // so then we should get adjacencies from mesh based geometry model.
370 : : //we should really check if (igeomIndex>=0);
371 [ + - ][ + - ]: 4156 : iGeom::Error err = igeom_instance()->getEntAdj(geom_handle(), (iBase_EntityType)dim, gents);
[ + - ]
372 [ + - ]: 4156 : IBERRCHK(err, "Trouble getting geom adjacencies.");
373 [ + - ]: 4156 : adjs.resize(gents.size());
374 : : iGeom::TagHandle mkmodeltag;
375 [ + - ][ + - ]: 4156 : err = igeom_instance()->getTagHandle("__MKModelEntityGeo", mkmodeltag);
376 [ + - ]: 4156 : IBERRCHK(err, "Failed to get tag handle for model entity.");
377 [ + - ][ + - ]: 8312 : err = igeom_instance()->getArrData(&gents[0], adjs.size(), mkmodeltag,
378 [ + - ][ + - ]: 8312 : &adjs[0]);
379 [ + - ]: 4156 : IBERRCHK(err, "Trouble getting ModelEnts for geom adjacencies.");
380 : 4156 : }
381 : :
382 : 94631 : int ModelEnt::dimension() const
383 : : {
384 [ + + ]: 94631 : if (iGeomEnt) {
385 : : iBase_EntityType tp;
386 [ + - ][ + - ]: 94587 : iGeom::Error err = igeom_instance()->getEntType(iGeomEnt, tp);
387 [ + - ]: 94587 : IBERRCHK(err, "Couldn't get geom entity type.");
388 : 94587 : return (tp - iBase_VERTEX);
389 : : }
390 [ + - ]: 44 : else if (moabEntSet) {
391 : : int dim;
392 [ + - ][ + - ]: 44 : moab::ErrorCode rval = moab_instance()->tag_get_data(mkCore->moab_geom_dim_tag(), &moabEntSet, 1, &dim);
[ + - ]
393 [ - + ][ # # ]: 44 : MBERRCHK(rval, moab_instance());
[ # # ][ # # ]
[ # # ]
394 : 44 : return dim;
395 : : }
396 [ # # ]: 94631 : else throw Error(MK_BAD_INPUT, "Couldn't get dimension for this ModelEnt.");
397 : : }
398 : :
399 : 672 : double ModelEnt::measure() const
400 : : {
401 : : double meas;
402 [ + - ][ + - ]: 672 : iGeom::Error err = igeom_instance()->measure(&iGeomEnt, 1, &meas);
403 [ + - ]: 672 : IBERRCHK(err, "Couldn't get geom entity measure.");
404 : 672 : return meas;
405 : : }
406 : :
407 : 0 : double ModelEnt::measure_discrete() const
408 : : {
409 : 0 : return measure();
410 : : }
411 : :
412 : 13488 : void ModelEnt::evaluate(double x, double y, double z,
413 : : double *close,
414 : : double *direction,
415 : : double *curvature1,
416 : : double *curvature2) const
417 : : {
418 : 13488 : iGeom::Error err = iBase_SUCCESS;
419 [ + + ][ - + ]: 13488 : if (0 == dimension() || 3 == dimension()) {
[ + + ]
420 [ + - ][ + - ]: 282 : if (direction || curvature1 || curvature2) {
[ - + ]
421 [ # # ][ # # ]: 0 : MKERRCHK(Error(MK_BAD_INPUT), "Direction or curvature not available for entities of this type.");
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
422 : : }
423 [ - + ]: 282 : else if (!close) {
424 [ # # ][ # # ]: 282 : MKERRCHK(Error(MK_BAD_INPUT), "Must pass closest point pointer for output.");
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
425 : : }
426 : : }
427 : :
428 [ + + ]: 13488 : if (0 == dimension()) {
429 : 282 : err = igeom_instance()->getVtxCoord(geom_handle(), close[0], close[1], close[2]);
430 : 282 : IBERRCHK(err, "Problem getting vertex coordinates.");
431 : 13488 : return;
432 : : }
433 [ - + ]: 13206 : else if (3 == dimension()) {
434 : 0 : err = igeom_instance()->getEntClosestPt(geom_handle(), x, y, z,
435 : 0 : close[0], close[1], close[2]);
436 : 0 : IBERRCHK(err, "Problem getting closest point for this model entity.");
437 : : }
438 : : else {
439 : : double cls[3], dir[3], curve1[3], curve2[3];
440 [ + + ]: 13206 : if (!close) close = cls;
441 [ + + ]: 13206 : if (!direction) direction = dir;
442 [ + - ]: 13206 : if (!curvature1) curvature1 = curve1;
443 [ + - ][ - + ]: 13206 : if (1 == dimension())
444 [ # # ]: 0 : err = igeom_instance()->getEgEvalXYZ(geom_handle(), x, y, z,
445 : 0 : close[0], close[1], close[2],
446 : 0 : direction[0], direction[1], direction[2],
447 [ # # ][ # # ]: 0 : curvature1[0], curvature1[1], curvature1[2]);
448 [ + - ][ + - ]: 13206 : else if (2 == dimension()) {
449 [ + - ]: 13206 : if (!curvature2) curvature2 = curve2;
450 [ + - ]: 13206 : err = igeom_instance()->getFcEvalXYZ(geom_handle(), x, y, z,
451 : 13206 : close[0], close[1], close[2],
452 : 13206 : direction[0], direction[1], direction[2],
453 : 13206 : curvature1[0], curvature1[1], curvature1[2],
454 [ + - ][ + - ]: 13206 : curvature2[0], curvature2[1], curvature2[2]);
455 : : }
456 : :
457 [ + - ]: 13206 : IBERRCHK(err, "Couldn't evaluate model entity.");
458 : : }
459 : : }
460 : :
461 : 997 : int ModelEnt::id() const
462 : : {
463 : : // get a global id for this entity, if there is one
464 [ + - ]: 997 : if (geom_handle()) {
465 : : iGeom::TagHandle gid_tag;
466 [ + - ][ + - ]: 997 : iGeom::Error err = igeom_instance()->getTagHandle("GLOBAL_ID", gid_tag);
467 [ + - ]: 997 : if (iBase_SUCCESS == err) {
468 : : int gid;
469 [ + - ][ + - ]: 997 : err = igeom_instance()->getIntData(geom_handle(), gid_tag, gid);
[ + - ]
470 [ + - ]: 997 : if (iBase_SUCCESS == err) return gid;
471 : : }
472 : : }
473 : :
474 [ # # ]: 0 : if (mesh_handle()) {
475 : : moab::Tag gid_tag;
476 : :
477 [ # # ][ # # ]: 0 : moab::ErrorCode err = mk_core()->moab_instance()->tag_get_handle("GLOBAL_ID",
478 : : 1, moab::MB_TYPE_INTEGER,
479 [ # # ]: 0 : gid_tag, moab::MB_TAG_SPARSE);
480 [ # # ]: 0 : if (moab::MB_SUCCESS == err) {
481 : : int gid;
482 [ # # ]: 0 : moab::EntityHandle this_mesh = mesh_handle();
483 [ # # ][ # # ]: 0 : err = mk_core()->moab_instance()->tag_get_data(gid_tag, &this_mesh, 1, &gid);
[ # # ]
484 [ # # ]: 0 : if (moab::MB_SUCCESS == err) return gid;
485 : : }
486 : : }
487 : :
488 : 997 : return -1;
489 : : }
490 : :
491 : 1938 : void ModelEnt::get_mesh(int dim,
492 : : std::vector<moab::EntityHandle> &ments,
493 : : bool bdy_too)
494 : : {
495 [ - + ]: 1938 : if (dim > dimension()) throw Error(MK_BAD_INPUT, "Called get_mesh for dimension %d on a %d-dimensional entity.",
496 [ # # ][ # # ]: 0 : dim, dimension());
497 : :
498 [ + + ][ + + ]: 1938 : if (!bdy_too || dim == dimension()) {
[ + + ]
499 : : // just owned entities, which will be in the set
500 : 1580 : moab::ErrorCode rval = mk_core()->moab_instance()->get_entities_by_dimension(moabEntSet, dimension(), ments);
501 [ - + ][ # # ]: 1580 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
502 : 1580 : return;
503 : : }
504 : :
505 : : // filter out the cases where order matters and boundary is requested too
506 [ + - ][ + - ]: 358 : else if (1 == dimension() && 0 == dim) {
[ + - ]
507 [ + - ][ + - ]: 716 : std::vector<moab::EntityHandle> tmp_edgevs, tmp_vvs;
508 [ + - ][ + - ]: 358 : moab::ErrorCode rval = mk_core()->moab_instance()->get_entities_by_dimension(moabEntSet, 0, tmp_edgevs);
[ + - ]
509 [ - + ][ # # ]: 358 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
510 [ + - ]: 358 : boundary(0, tmp_vvs);
511 [ - + ][ # # ]: 358 : if (tmp_vvs.empty()) throw Error(MK_NOT_FOUND, "No mesh on bounding entity.");
512 [ + - ][ + - ]: 358 : ments.push_back(tmp_vvs[0]);
513 [ + - ][ + - ]: 358 : std::copy(tmp_edgevs.begin(), tmp_edgevs.end(), std::back_inserter(ments));
514 [ + + ][ + - ]: 358 : if (2 == tmp_vvs.size()) ments.push_back(tmp_vvs[1]);
[ + - ]
515 [ + - ][ + - ]: 24 : else if (1 == tmp_vvs.size()) ments.push_back(tmp_vvs[0]);
[ + - ]
516 : 716 : return;
517 : : }
518 : : else {
519 : : // remaining case is bdy_too=true and dim < dimension()
520 : : // first, get dimension()-dimensional entities, then get dim-dimensional neighbors of those
521 [ # # ]: 0 : std::vector<moab::EntityHandle> tmp_ments;
522 [ # # ][ # # ]: 0 : moab::ErrorCode rval = mk_core()->moab_instance()->get_entities_by_dimension(moabEntSet, dimension(), tmp_ments);
[ # # ][ # # ]
523 [ # # ][ # # ]: 0 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
524 [ # # ][ # # ]: 0 : rval = mk_core()->moab_instance()->get_adjacencies(&tmp_ments[0], tmp_ments.size(), dim, false, ments,
[ # # ]
525 [ # # ]: 0 : moab::Interface::UNION);
526 [ # # ][ # # ]: 1938 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
527 : : }
528 : : }
529 : :
530 : 740 : void ModelEnt::boundary(int dim,
531 : : std::vector<moab::EntityHandle> &bdy,
532 : : std::vector<int> *senses,
533 : : std::vector<int> *group_sizes)
534 : : {
535 [ + - ]: 740 : MEntVector me_loop;
536 [ + - ][ + - ]: 1480 : std::vector<int> me_senses, me_group_sizes;
537 [ + - ]: 740 : boundary(dim, me_loop, &me_senses, &me_group_sizes);
538 : :
539 : 740 : MEntVector::iterator vit = me_loop.begin();
540 : 740 : std::vector<int>::iterator sense_it = me_senses.begin(), grpsize_it = me_group_sizes.begin();
541 : 740 : int gents_ctr = 0, first_ment = 0;
542 : :
543 : : // packing into a single list, but need to keep track of loop/shell sizes; don't increment
544 : : // grpsize_it here, just when we've done all the mes in this group
545 [ + - ][ + - ]: 2472 : for (vit = me_loop.begin(); vit != me_loop.end(); vit++, sense_it++) {
[ + - ][ + + ]
546 [ + - ]: 1732 : std::vector<moab::EntityHandle> tmp_ments;
547 [ + - ][ + - ]: 1732 : (*vit)->get_mesh(dim, tmp_ments, true);
548 [ + - ][ + + ]: 1732 : if (*sense_it == SENSE_REVERSE)
549 [ + - ]: 197 : std::reverse(tmp_ments.begin(), tmp_ments.end());
550 [ + - ][ + + ]: 1732 : if (2 == dimension() && 0 == dim)
[ + + ][ + + ]
551 : : // assembling vertices on loops; don't do last, since that'll be the first of the
552 : : // next edge
553 [ + - ][ + - ]: 284 : std::copy(tmp_ments.begin(), tmp_ments.end()-1, std::back_inserter(bdy));
[ + - ]
554 : : else
555 [ + - ][ + - ]: 1448 : std::copy(tmp_ments.begin(), tmp_ments.end(), std::back_inserter(bdy));
556 [ + + ]: 1732 : if (senses) {
557 [ + - ][ + - ]: 12546 : for (unsigned int i = 0; i < tmp_ments.size(); i++) senses->push_back(*sense_it);
[ + + ]
558 : : }
559 [ + + ]: 1732 : if (group_sizes) {
560 : 328 : gents_ctr++;
561 [ + - ][ + + ]: 328 : if (gents_ctr == *grpsize_it) {
562 [ + - ]: 92 : group_sizes->push_back(bdy.size()-first_ment);
563 : : // first_ment will be the first element in the next group
564 : 92 : first_ment = bdy.size();
565 : 92 : gents_ctr = 0;
566 [ + - ]: 92 : grpsize_it++;
567 : : } // block at end of a given loop
568 : : } // block updating group size
569 : 2472 : } // over loop
570 : :
571 : 740 : }
572 : :
573 : 820 : void ModelEnt::boundary(int dim,
574 : : MEntVector &entities,
575 : : std::vector<int> *senses,
576 : : std::vector<int> *group_sizes)
577 : : {
578 : : // for a given entity, return the bounding edges in the form of edge groups,
579 : : // oriented ccw around surface
580 [ + - ]: 820 : int this_dim = dimension();
581 : :
582 : : // shouldn't be calling this function if we're a vertex
583 [ - + ][ # # ]: 820 : if (this_dim == iBase_VERTEX) throw Error(MK_BAD_INPUT, "Shouldn't call this for a vertex.");
584 [ - + ][ # # ]: 820 : else if (dim >= this_dim) throw Error(MK_BAD_INPUT, "Calling for boundary entities of equal or greater dimension.");;
585 : :
586 : : // get (d-1)-adj ents & senses
587 [ + - ]: 820 : MEntVector tmp_ents;
588 [ + - ]: 820 : get_adjacencies(this_dim-1, tmp_ents);
589 : :
590 : : // get out here if we're a curve
591 [ + + ]: 820 : if (1 == this_dim) {
592 : : // it is important to order the vertices, especially for boundary
593 : : // cases get_adjacency does not say anything about order
594 [ + + ]: 645 : if (tmp_ents.size() <= 1) // no worry
595 : : {
596 [ + + ]: 80 : for (unsigned int i = 0; i < tmp_ents.size(); i++) {
597 [ + - ][ + - ]: 40 : entities.push_back(tmp_ents[i]);
598 [ + - ][ + - ]: 40 : if (senses) senses->push_back(SENSE_FORWARD);
599 : : }
600 : : }
601 : : else
602 : : {
603 : : // we have at least 2 vertices; if more than 2, this is an error
604 [ - + ]: 605 : if (tmp_ents.size() > 2)
605 [ # # ]: 0 : throw Error(MK_FAILURE, " edge with too many vertices ");
606 : : // we have now exactly 2 vertices; order them according to the
607 : : // edge they come from
608 : 605 : int sense = 0;
609 [ + - ][ + - ]: 1210 : iGeom::Error rg = igeom_instance()->getEgVtxSense( geom_handle(),tmp_ents[0]->geom_handle(),
610 [ + - ][ + - ]: 1210 : tmp_ents[1]->geom_handle(), sense );
[ + - ][ + - ]
[ + - ]
611 [ + - ]: 605 : IBERRCHK(rg, "Trouble getting edge sense");
612 [ + + ]: 605 : if (-1==sense)
613 : : {
614 : : // the vertices are reversed in adjacency list
615 [ + - ][ + - ]: 25 : entities.push_back(tmp_ents[1]);
616 [ + - ][ + - ]: 25 : entities.push_back(tmp_ents[0]);
617 : : }
618 : : else
619 : : {
620 : : // vertices are fine
621 [ + - ][ + - ]: 580 : entities.push_back(tmp_ents[0]);
622 [ + - ][ + - ]: 580 : entities.push_back(tmp_ents[1]);
623 : : }
624 [ + - ]: 605 : if (senses) {
625 [ + - ]: 605 : senses->push_back(SENSE_FORWARD);
626 [ + - ]: 605 : senses->push_back(SENSE_FORWARD);
627 : : }
628 : : }
629 : 645 : return; // we need to get out, we treated a dim 1 entity
630 : : }
631 : :
632 : : // get adjacent entities into a sorted, mutable list
633 [ + - ][ + + ]: 995 : MEntVector b_ents = tmp_ents;
634 [ + - ][ + - ]: 350 : MEntSet dbl_curves, sgl_curves;
[ + - ][ + - ]
635 : :
636 [ + - ]: 175 : std::sort(b_ents.begin(), b_ents.end());
637 : :
638 [ + + ][ + - ]: 392 : while (!b_ents.empty()) {
639 : : // get 1st in group
640 : :
641 [ + - ][ + - ]: 434 : MEntVector group_stack, this_group;
[ + - ]
642 [ + - ][ + - ]: 217 : group_stack.push_back(b_ents.front());
643 : :
644 : 217 : int current_senses_size = 0;
645 [ + - ]: 217 : if (senses)
646 : 217 : current_senses_size= senses->size();
647 : : // while there are still entities on the stack
648 [ + + ]: 1091 : while (!group_stack.empty()) {
649 : :
650 : : // pop one off & get its sense
651 [ + - ][ + - ]: 874 : ModelEnt *this_entity = group_stack.back(); group_stack.pop_back();
652 : :
653 : : int this_sense;
654 : : iGeom::Error err =
655 [ + - ][ + - ]: 874 : igeom_instance()->getSense(this_entity->geom_handle(), geom_handle(), this_sense);
[ + - ][ + - ]
656 [ + - ]: 874 : IBERRCHK(err,"Error getting geometry sense");
657 : :
658 : : // if we already have this one, continue
659 [ + + ][ + + ]: 2582 : if ((0 != this_sense &&
660 [ + + ][ + - ]: 3476 : std::find(this_group.begin(), this_group.end(), this_entity) != this_group.end()) ||
[ + - ][ - + ]
[ + + ][ + + ]
[ # # # # ]
661 [ + - ][ + - ]: 1676 : std::find(dbl_curves.begin(), dbl_curves.end(), this_entity) != dbl_curves.end())
[ + + ][ + + ]
[ # # # # ]
662 : 72 : continue;
663 : :
664 : : // either way we need the d-2 entities
665 [ + - ]: 802 : MEntVector bridges;
666 [ + - ][ + - ]: 802 : this_entity->get_adjacencies(dimension()-2, bridges);
667 [ + - ][ + + ]: 802 : if (dimension()-2==0 && bridges.size()==2)// edge case; this_entity is an edge
[ + + ][ + + ]
668 : : {
669 : : // if the vertices are reversed, change them back
670 : 710 : int sense_vertices = 0;
671 [ + - ]: 710 : iGeom::Error rg = igeom_instance()->getEgVtxSense(this_entity->geom_handle(),
672 [ + - ]: 710 : bridges[0]->geom_handle(),
673 [ + - ][ + - ]: 1420 : bridges[1]->geom_handle(), sense_vertices );
[ + - ][ + - ]
[ + - ]
674 [ + - ]: 710 : IBERRCHK(rg, "Trouble getting edge sense");
675 [ + + ]: 710 : if (-1==sense_vertices)
676 : : {
677 : : // reverse the order of bridges;
678 : : // a better fix would be to have the bridges from get_adjacencies in order
679 : : // this would mean that children of sets should be returned in the order they
680 : : // were inserted as children (maybe too hard to control that)
681 [ + - ]: 18 : ModelEnt * tmp = bridges[0];
682 [ + - ][ + - ]: 18 : bridges[0] = bridges[1];
683 [ + - ]: 710 : bridges[1] = tmp;
684 : : }
685 : :
686 : : }
687 : :
688 : : // only remove from the list of candidates if it's not dual-sensed
689 [ + + ]: 802 : if (0 != this_sense)
690 [ + - ][ + - ]: 782 : b_ents.erase(std::remove(b_ents.begin(), b_ents.end(), this_entity), b_ents.end());
691 : :
692 : : // if it's double-sensed and this is the first time we're seeing it, find the right sense
693 : : else {
694 [ + - ][ - + ]: 20 : assert(dimension() == 2);
695 [ + - ][ + - ]: 20 : if (std::find(sgl_curves.begin(), sgl_curves.end(), this_entity) == sgl_curves.end()) {
[ + + ]
696 [ + - ]: 10 : sgl_curves.insert(this_entity);
697 [ + + ]: 10 : if (!this_group.empty()) {
698 [ + - ][ + - ]: 5 : ModelEnt *common_v = this_entity->shared_entity(this_group.back(), 0);
699 [ + - ][ - + ]: 5 : if (common_v == bridges[0]) this_sense = 1;
700 [ + - ][ + - ]: 5 : else if (common_v == bridges[1]) this_sense = -1;
701 : 5 : else assert(false);
702 : : }
703 : : // else, if this is the first one in the loop, just choose a sense
704 : : else {
705 : 10 : this_sense = 1;
706 : : }
707 : : }
708 : : else {
709 : : // else this is the second time we're seeing it, move it to the dbl_curves list and remove
710 : : // from the candidates list
711 [ + - ]: 10 : dbl_curves.insert(this_entity);
712 [ + - ]: 10 : sgl_curves.erase(this_entity);
713 [ + - ][ + - ]: 10 : b_ents.erase(std::remove(b_ents.begin(), b_ents.end(), this_entity), b_ents.end());
714 [ + - ]: 10 : if (senses) {
715 : : // need to get the sense of the first appearance in list
716 [ + - ]: 10 : MEntVector::iterator vit = std::find(this_group.begin(), this_group.end(), this_entity);
717 [ + - ][ - + ]: 10 : assert(vit != this_group.end());
718 [ + - ][ + - ]: 10 : int other_sense = (*senses)[current_senses_size + vit-this_group.begin()];
[ + - ]
719 : 10 : this_sense = SENSE_REVERSE*other_sense;
720 : : }
721 : : }
722 : : }
723 : :
724 : : // it's in the group; put on group & remove from untreated ones
725 [ + - ]: 802 : this_group.push_back(this_entity);
726 [ + - ][ + - ]: 802 : if (senses) senses->push_back(this_sense);
727 [ + - ][ + - ]: 1604 : MEntVector tmp_from, tmp_adjs;
[ + - ][ + - ]
728 : :
729 : : // if we're on a face and we're the first in a group, check sense of this first
730 : : // edge; make sure "next" in loop sense is last on list
731 [ + - ][ + + ]: 802 : if (2 == dimension() && this_group.size() == 1) {
[ + + ][ + + ]
732 : :
733 : : // get vertex which we know is shared by the "right" next edge; first get the vertices
734 : : // if sense of current edge is forward, it's the 2nd vertex we want,
735 : : // otherwise the first
736 [ + + ][ + + ]: 213 : if (1 == this_sense && bridges.size() > 1) tmp_from.push_back(bridges[1]);
[ + + ][ + - ]
[ + - ]
737 [ + - ][ + - ]: 91 : else tmp_from.push_back(bridges[0]);
738 [ + - ]: 213 : tmp_adjs = b_ents;
739 [ + - ][ + - ]: 213 : get_adjs_bool(tmp_from, dimension()-1, tmp_adjs, INTERSECT);
740 : : }
741 : : else {
742 [ + - ][ + - ]: 589 : std::copy(bridges.begin(), bridges.end(), std::back_inserter(tmp_from));
743 [ + - ]: 589 : MEntVector tmp_adjs2;
744 [ + - ][ + - ]: 589 : get_adjs_bool(tmp_from, dimension()-1, tmp_adjs2, UNION);
745 [ + - ]: 589 : std::sort(tmp_adjs2.begin(), tmp_adjs2.end());
746 [ + - ]: 589 : tmp_adjs.resize(tmp_adjs2.size());
747 : : tmp_adjs.erase(std::set_intersection(b_ents.begin(), b_ents.end(),
748 : : tmp_adjs2.begin(), tmp_adjs2.end(),
749 [ + - ][ + - ]: 589 : tmp_adjs.begin()), tmp_adjs.end());
750 : : }
751 : :
752 [ + - ][ + + ]: 802 : if (2 == dimension() && tmp_adjs.size() > 1) {
[ + + ][ + + ]
753 : : // more than one adjacent edge - need to evaluate winding angle to find right one
754 [ + - ]: 10 : ModelEnt *next_ent = next_winding(this_entity, this_sense, tmp_adjs);
755 [ - + ]: 10 : if (NULL == next_ent) return;
756 [ + - ]: 10 : group_stack.push_back(next_ent);
757 : : }
758 [ + + ]: 792 : else if (!tmp_adjs.empty())
759 [ + - ][ + - ]: 802 : std::copy(tmp_adjs.begin(), tmp_adjs.end(), std::back_inserter(group_stack));
[ + - ]
760 : 802 : }
761 : :
762 : : // put group in group list
763 [ + - ][ + - ]: 217 : std::copy(this_group.begin(), this_group.end(), std::back_inserter(entities));
764 [ + + ][ + - ]: 217 : if (group_sizes) group_sizes->push_back(this_group.size());
[ + - ]
765 : 1037 : }
766 : : }
767 : :
768 : 5 : ModelEnt *ModelEnt::shared_entity(ModelEnt *ent2, int to_dim)
769 : : {
770 : : // find the shared entity between the two entities of the prescribed dimension
771 [ + - ][ + - ]: 10 : MEntVector from_ents(2), to_ents;
772 [ + - ]: 5 : from_ents[0] = this;
773 [ + - ]: 5 : from_ents[1] = ent2;
774 : : try {
775 [ + - ]: 5 : get_adjs_bool(from_ents, to_dim, to_ents, INTERSECT, false);
776 : : }
777 [ # # # # ]: 0 : catch (Error err) {
778 [ # # # # : 0 : if (err.error_code() != MK_SUCCESS || to_ents.empty()) return NULL;
# # # # #
# ]
779 : 0 : }
780 : :
781 [ - + ][ # # ]: 5 : if (to_ents.size() > 1) throw Error(MK_MULTIPLE_FOUND, "Multiple shared entities found.");
782 : :
783 [ + - ]: 10 : return to_ents[0];
784 : : }
785 : :
786 : 807 : void ModelEnt::get_adjs_bool(MEntVector &from_ents,
787 : : int to_dim,
788 : : MEntVector &to_ents,
789 : : BooleanType op_type,
790 : : bool only_to_ents)
791 : : {
792 [ - + ]: 807 : if (from_ents.empty()) {
793 : 0 : to_ents.clear();
794 : 807 : return;
795 : : }
796 : :
797 [ + - ]: 807 : MEntVector bridges;
798 : :
799 : 807 : MEntVector::iterator from_it = from_ents.begin();
800 [ + + ][ + + ]: 807 : if (to_ents.empty() && !only_to_ents && op_type == INTERSECT) {
[ + - ][ + + ]
801 [ + - ][ + - ]: 5 : from_ents.front()->get_adjacencies(to_dim, bridges);
802 [ + - ][ + - ]: 5 : std::copy(bridges.begin(), bridges.end(), std::back_inserter(to_ents));
803 [ + - ]: 5 : from_it++;
804 : : }
805 : :
806 [ + - ]: 807 : std::sort(to_ents.begin(), to_ents.end());
807 [ + - ]: 1614 : MEntVector result_ents(to_ents.size());
808 [ + - ][ + - ]: 2348 : for (; from_it != from_ents.end(); from_it++) {
[ + + ]
809 : 1541 : bridges.clear();
810 [ + - ][ + - ]: 1541 : (*from_it)->get_adjacencies(to_dim, bridges);
811 [ + + ]: 1541 : if (op_type == INTERSECT) {
812 [ + - ]: 218 : std::sort(bridges.begin(), bridges.end());
813 : : result_ents.erase(std::set_intersection(to_ents.begin(), to_ents.end(),
814 : : bridges.begin(), bridges.end(),
815 [ + - ][ + - ]: 218 : result_ents.begin()), result_ents.end());
816 : : }
817 : : else {
818 [ + - ][ + - ]: 1323 : std::copy(bridges.begin(), bridges.end(), std::back_inserter(result_ents));
819 : : }
820 : :
821 [ + - ]: 1541 : to_ents = result_ents;
822 : : }
823 : :
824 [ + + ]: 807 : if (op_type == UNION)
825 [ + - ]: 589 : std::sort(to_ents.begin(), to_ents.end());
826 : :
827 [ + - ][ + - ]: 1614 : to_ents.erase(std::unique(to_ents.begin(), to_ents.end()), to_ents.end());
828 : : }
829 : :
830 : 10 : ModelEnt *ModelEnt::next_winding(ModelEnt *this_edge,
831 : : int this_sense,
832 : : MEntVector &tmp_adjs)
833 : : {
834 : : // given this_entity, a d-1 entity bounding gentity, and optional "next"
835 : : // entities in tmp_adjs, find the next one based on windings around the shared
836 : : // vertex
837 : :
838 : : // first, get the shared vertex
839 [ + - ]: 10 : MEntVector verts;
840 [ + - ]: 10 : this_edge->get_adjacencies(0, verts);
841 : :
842 [ + - ]: 10 : ModelEnt *shared_vert = verts[0];
843 [ + + ][ + - ]: 10 : if (this_sense == 1 && verts.size() > 1) shared_vert = verts[1];
[ + + ][ + - ]
844 : :
845 : : // get locations just before the vertex, at the vertex, and just after the vertex
846 : : double v1[3], v2[3], v3[3];
847 : : double umin, umax;
848 [ + - ][ + - ]: 10 : igeom_instance()->getEntURange(this_edge->geom_handle(), umin, umax);
[ + - ]
849 : : double utgt;
850 [ + + ]: 10 : if (1 == this_sense) utgt = umin + 0.9 * (umax - umin);
851 : 5 : else utgt = umin + 0.1 * (umax - umin);
852 [ + - ][ + - ]: 10 : igeom_instance()->getEntUtoXYZ(this_edge->geom_handle(), utgt, v1[0], v1[1], v1[2]);
[ + - ]
853 [ + - ][ + - ]: 10 : igeom_instance()->getVtxCoord(shared_vert->geom_handle(), v2[0], v2[1], v2[2]);
[ + - ]
854 : 10 : double v21[] = {v1[0]-v2[0], v1[1]-v2[1], v1[2]-v2[2]};
855 [ + - ]: 10 : VecUtil::normalize(v21);
856 : :
857 : : // get the normal vector at the vertex
858 : : double normal[3];
859 [ + - ]: 10 : igeom_instance()->getEntNrmlXYZ(geom_handle(), v2[0], v2[1], v2[2],
860 [ + - ][ + - ]: 10 : normal[0], normal[1], normal[2]);
861 : :
862 : : // now loop over candidates, finding magnitude of swept angle
863 : 10 : ModelEnt *other = NULL;
864 : 10 : double angle = VecUtil::TWO_PI;
865 : : bool is_adj;
866 [ + - ][ + - ]: 30 : for (MEntVector::iterator vit = tmp_adjs.begin(); vit != tmp_adjs.end(); vit++) {
[ + + ]
867 : : // if we're here, we have multiple candidates, therefore don't choose the same one
868 [ + - ][ + + ]: 20 : if (*vit == this_edge)
869 : 10 : continue;
870 [ + - ][ + - ]: 10 : igeom_instance()->isEntAdj((*vit)->geom_handle(), shared_vert->geom_handle(), is_adj);
[ + - ][ + - ]
[ + - ]
871 [ - + ]: 10 : if (!is_adj)
872 : 0 : continue;
873 : :
874 : : // get param range
875 [ + - ][ + - ]: 10 : igeom_instance()->getEntURange((*vit)->geom_handle(), umin, umax);
[ + - ][ + - ]
876 : : // get sense
877 : : int tmp_sense;
878 [ + - ][ + - ]: 10 : igeom_instance()->getSense((*vit)->geom_handle(), geom_handle(), tmp_sense);
[ + - ][ + - ]
[ + - ]
879 [ + + ]: 10 : if (1 == tmp_sense) utgt = umin + 0.1 * (umax - umin);
880 : 5 : else utgt = umin + 0.9 * (umax - umin);
881 [ + - ][ + - ]: 10 : igeom_instance()->getEntUtoXYZ((*vit)->geom_handle(), utgt, v3[0], v3[1], v3[2]);
[ + - ][ + - ]
882 : 10 : double v23[] = {v3[0]-v2[0], v3[1]-v2[1], v3[2]-v2[2]};
883 [ + - ]: 10 : VecUtil::normalize(v23);
884 : :
885 [ + - ]: 10 : VecUtil::cross(normal, v23, v1);
886 : :
887 [ + - ]: 10 : VecUtil::cross(v1, normal, v3);
888 : :
889 [ + - ]: 10 : double x = VecUtil::dot(v21, v3);
890 [ + - ]: 10 : double y = VecUtil::dot(v21, v1);
891 : :
892 [ + + ][ - + ]: 10 : assert(x != 0.0 || y != 0.0);
893 : 10 : double this_angle = atan2(y, x);
894 [ - + ]: 10 : if (this_angle < 0.0)
895 : : {
896 : 0 : this_angle += VecUtil::TWO_PI;
897 : : }
898 [ + - ]: 10 : if (this_angle < angle) {
899 [ + - ]: 10 : other = *vit;
900 : 10 : angle = this_angle;
901 : : }
902 : : }
903 : :
904 : 10 : return other;
905 : : }
906 : :
907 : 0 : void ModelEnt::boundary(int dim,
908 : : moab::Range &ents) const
909 : : {
910 [ # # ]: 0 : throw Error(MK_NOT_IMPLEMENTED, "Not implemented yet.");
911 : : }
912 : :
913 : 76 : void ModelEnt::get_indexed_connect_coords(std::vector<moab::EntityHandle> &ents,
914 : : std::vector<int> *senses,
915 : : moab::Tag tagh,
916 : : std::vector<int> &ents_as_ids,
917 : : std::vector<double> &coords,
918 : : moab::Range *verts_range,
919 : : int start_index)
920 : : {
921 : : // if we need to, make a tag
922 : : moab::ErrorCode rval;
923 : 76 : bool i_created_tag = false;
924 [ + - ]: 76 : if (0 == tagh) {
925 : : //rval = mk_core()->moab_instance()->tag_create("__ModelEntidtag", sizeof(int), moab::MB_TAG_DENSE,
926 : : // moab::MB_TYPE_INTEGER, tagh, NULL);
927 [ + - ][ + - ]: 76 : rval = mk_core()->moab_instance()->tag_get_handle("__ModelEntidtag", 1, moab::MB_TYPE_INTEGER, tagh,
928 [ + - ]: 76 : moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);
929 [ - + ][ # # ]: 76 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
930 : 76 : i_created_tag = true;
931 : : }
932 : :
933 : : // put entities into range, after clearing it
934 [ + - ][ + - ]: 152 : moab::Range tmp_range, ents_range;
935 [ - + ]: 76 : if (!verts_range) verts_range = &tmp_range;
936 [ + - ]: 76 : verts_range->clear();
937 : :
938 [ + - ][ + - ]: 76 : std::copy(ents.begin(), ents.end(), moab::range_inserter(ents_range));
939 [ + - ][ + - ]: 228 : bool all_verts = (mkCore->moab_instance()->type_from_handle(*ents_range.begin()) ==
[ + - ][ + - ]
[ - + ][ # # ]
940 [ + - ][ + - ]: 228 : mkCore->moab_instance()->type_from_handle(*ents_range.begin()) &&
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ][ # # ]
941 [ + - ][ + - ]: 228 : mkCore->moab_instance()->type_from_handle(*ents_range.begin()) == moab::MBVERTEX);
[ + - ][ + - ]
[ + - ][ # # ]
942 : :
943 : :
944 : : // get connectivity of all ents and store in range
945 [ + + ]: 76 : if (!all_verts) {
946 [ + - ][ + - ]: 34 : rval = mkCore->moab_instance()->get_connectivity(ents_range, *verts_range);
947 [ - + ][ # # ]: 34 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
948 : : }
949 [ + - ]: 42 : else *verts_range = ents_range;
950 : :
951 : : // resize index array, to max number of connectivity entries
952 [ + - ][ + - ]: 76 : int max_numconnect = moab::CN::VerticesPerEntity(mkCore->moab_instance()->type_from_handle(*(ents_range.rbegin())));
[ + - ][ + - ]
[ + - ]
953 [ + - ]: 76 : ents_as_ids.resize(ents.size()*max_numconnect);
954 : :
955 : : // temporarily store ids in ents_as_ids array, and set id tag
956 [ + - ][ - + ]: 76 : assert(ents_as_ids.size() >= verts_range->size());
957 [ + - ][ + - ]: 8674 : for (unsigned int i = 0; i < verts_range->size(); i++) ents_as_ids[i] = i+start_index;
[ + + ]
958 [ + - ][ + - ]: 76 : rval = mk_core()->moab_instance()->tag_set_data(tagh, *verts_range, &ents_as_ids[0]);
[ + - ][ + - ]
959 [ - + ][ # # ]: 76 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
960 : :
961 : : // get ids into ids vector in connectivity or ents order
962 [ + + ]: 76 : if (all_verts) {
963 [ + - ][ + - ]: 42 : rval = mk_core()->moab_instance()->tag_get_data(tagh, &ents[0], ents.size(), &ents_as_ids[0]);
[ + - ][ + - ]
[ + - ]
964 [ - + ][ # # ]: 42 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
965 [ + - ]: 42 : ents_as_ids.resize(ents.size());
966 : : }
967 : : else {
968 : : // loop over entities
969 [ + - ]: 34 : std::vector<moab::EntityHandle> storage;
970 : : const moab::EntityHandle *connect;
971 : : int num_connect;
972 : 34 : unsigned int last = 0;
973 [ + + ]: 12298 : for (unsigned int i = 0; i < ents.size(); i++) {
974 : : // first, get connect vector
975 [ + - ][ + - ]: 12264 : mk_core()->moab_instance()->get_connectivity(ents[i], connect, num_connect, true, &storage);
[ + - ][ + - ]
976 [ - + ][ # # ]: 12264 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
977 [ - + ]: 12264 : assert(ents_as_ids.size() >= last + num_connect);
978 : : // next, get ids and put into ids vector
979 [ + - ][ + - ]: 12264 : rval = mk_core()->moab_instance()->tag_get_data(tagh, connect, num_connect, &ents_as_ids[last]);
[ + - ][ + - ]
980 [ - + ][ # # ]: 12264 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
981 : : // reverse if necessary
982 [ + - ][ + - ]: 12264 : assert(!senses || (*senses)[i] != SENSE_BOTH);
[ - + ]
983 [ + - ][ + - ]: 12264 : if (senses && (*senses)[i] == SENSE_REVERSE) std::reverse(&ents_as_ids[last], &ents_as_ids[last+num_connect]);
[ + + ][ + + ]
[ + - ][ + - ]
[ + - ]
984 : : // update last
985 : 12264 : last += num_connect;
986 : 34 : }
987 : : }
988 : :
989 : : // get coords for range-ordered vertices
990 [ + - ][ + - ]: 76 : coords.resize(3*verts_range->size());
991 [ + - ][ + - ]: 76 : rval = mk_core()->moab_instance()->get_coords(*verts_range, &coords[0]);
[ + - ][ + - ]
992 [ - + ][ # # ]: 76 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
993 : :
994 : : // delete the tag if I created it
995 [ + - ]: 76 : if (i_created_tag) {
996 [ + - ][ + - ]: 76 : rval = mk_core()->moab_instance()->tag_delete(tagh);
[ + - ]
997 [ - + ][ # # ]: 76 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
998 : 76 : }
999 : 76 : }
1000 : :
1001 : 0 : iGeom::EntityHandle ModelEnt::geom_handle(moab::EntityHandle ment) const
1002 : : {
1003 : : // use iRel to get this information or use model entity tag
1004 : 0 : iBase_EntityHandle gent = 0;
1005 [ # # ]: 0 : if (-1 != irelIndex)
1006 : : {
1007 [ # # ][ # # ]: 0 : Error err = mkCore->irel_pair(irelIndex)->getSetEntRelation(IBSH(ment), true, gent);
[ # # ]
1008 [ # # ][ # # ]: 0 : MKERRCHK(err, "Failed to get geometry handle for mesh set.");
[ # # ][ # # ]
[ # # ]
1009 : 0 : return gent;
1010 : : }
1011 : : else
1012 : : {
1013 : : // use the model entity tag from moab instance... assume only one here
1014 : : //moab::Tag moabModelTag = mkCore->moab_model_tag();
1015 : 0 : ModelEnt *modelEnt = NULL;
1016 [ # # ][ # # ]: 0 : moab::ErrorCode rval = mkCore->moab_instance()->tag_get_data(mkCore->moab_model_tag(), &ment, 1, &modelEnt);
[ # # ]
1017 [ # # ][ # # ]: 0 : MBERRCHK(rval, mkCore->moab_instance());
[ # # ][ # # ]
[ # # ]
1018 [ # # ]: 0 : if (NULL == modelEnt)
1019 : 0 : return gent;// still null
1020 [ # # ]: 0 : return modelEnt->geom_handle();
1021 : : }
1022 : : }
1023 : :
1024 : 135999 : iGeom *ModelEnt::igeom_instance() const
1025 : : {
1026 [ - + ][ # # ]: 135999 : if (-1 == igeomIndex) throw Error(MK_FAILURE, "No geometry instance associated with this ModelEnt.");
1027 : 135999 : return mkCore->igeom_instance(igeomIndex);
1028 : : }
1029 : :
1030 : 44 : moab::Interface *ModelEnt::moab_instance() const
1031 : : {
1032 [ - + ][ # # ]: 44 : if (-1 == meshIndex) throw Error(MK_FAILURE, "No moab instance associated with this ModelEnt.");
1033 : 44 : return mkCore->moab_instance(meshIndex);
1034 : : }
1035 : :
1036 : 3722 : moab::EntityHandle ModelEnt::mesh_handle(iGeom::EntityHandle gent) const
1037 : : {
1038 : : // use iRel to get this information or use model entity tag
1039 : 3722 : moab::EntityHandle ment = 0;
1040 [ + + ]: 3722 : if (-1 != irelIndex)
1041 : : {
1042 : : iBase_EntitySetHandle h;
1043 [ + - ][ + - ]: 3590 : iRel::Error err = mkCore->irel_pair(irelIndex)->getEntSetRelation(gent, false, h);
1044 [ + + ]: 3590 : IBERRCHK(err, "Failed to get mesh set handle for geometry entity.");
1045 : 1795 : return reinterpret_cast<moab::EntityHandle>(h);
1046 : : }
1047 : : else
1048 : : {
1049 : : // if no irel, use the model entity tag from the igeom instance
1050 : : //
1051 [ + - ][ - + ]: 132 : if (!igeom_instance())
1052 : 0 : return ment;// NULL so far
1053 : : iGeom::TagHandle mkmodeltag;
1054 [ + - ][ + - ]: 132 : iBase_ErrorType err = igeom_instance()->getTagHandle("__MKModelEntityGeo", mkmodeltag);
1055 [ + - ]: 132 : IBERRCHK(err, "Failed to get tag handle for model entity.");
1056 : : // now get the model entity
1057 : 132 : ModelEnt * modelEnt = NULL;
1058 [ + - ][ + - ]: 132 : err = igeom_instance()->getData(gent, mkmodeltag, &modelEnt);
1059 [ + + ][ - + ]: 132 : if (NULL == modelEnt || iBase_TAG_NOT_FOUND == err)
1060 : 66 : return ment;// still null
1061 [ + - ]: 3722 : return modelEnt->mesh_handle();
1062 : : }
1063 : :
1064 : : }
1065 : :
1066 : 0 : moab::EntityHandle ModelEnt::mesh_handle(iGeom::EntitySetHandle gent) const
1067 : : {
1068 : : // use iRel to get this information or use model entity tag
1069 : 0 : moab::EntityHandle ment = 0;
1070 [ # # ]: 0 : if (-1 != irelIndex)
1071 : : {
1072 : : iBase_EntitySetHandle h;
1073 [ # # ][ # # ]: 0 : iRel::Error err = mkCore->group_set_pair(irelIndex)->getSetSetRelation(gent, false, h);
1074 [ # # ]: 0 : IBERRCHK(err, "Failed to get mesh set handle for geometry entity.");
1075 : 0 : return reinterpret_cast<moab::EntityHandle>(h);
1076 : : }
1077 : : else
1078 : : {
1079 : : // if no irel, use the model entity tag from the igeom instance
1080 : : //
1081 [ # # ]: 0 : if (-1 == igeomIndex)
1082 : 0 : return ment;// NULL so far
1083 : : iGeom::TagHandle mkmodeltag;
1084 [ # # ][ # # ]: 0 : iBase_ErrorType err = igeom_instance()->getTagHandle("__MKModelEntityGeo", mkmodeltag);
1085 [ # # ]: 0 : IBERRCHK(err, "Failed to get tag handle for model entity.");
1086 : : // now get the model entity
1087 : 0 : ModelEnt * modelEnt = NULL;
1088 [ # # ][ # # ]: 0 : err = igeom_instance()->getEntSetData(gent, mkmodeltag, &modelEnt);
1089 [ # # ][ # # ]: 0 : if (NULL == modelEnt || iBase_TAG_NOT_FOUND == err)
1090 : 0 : return ment;// still null
1091 [ # # ]: 0 : return modelEnt->mesh_handle();
1092 : : }
1093 : :
1094 : : }
1095 : :
1096 : : /** \brief Get mesh interval size, if any
1097 : : * Returns -1 if no size set on this entity. If intervals are set, returns computed size.
1098 : : * \return Interval size for this ModelEnt.
1099 : : */
1100 : 44 : double ModelEnt::mesh_interval_size() const
1101 : : {
1102 [ + - ]: 88 : if (-1 != sizing_function_index() && mk_core()->sizing_function(sizing_function_index()) &&
[ + - + - ]
[ + - ]
1103 : 44 : -1 != mk_core()->sizing_function(sizing_function_index())->size())
1104 : 44 : return mk_core()->sizing_function(sizing_function_index())->size();
1105 : :
1106 [ # # ][ # # ]: 0 : else if (1 == dimension() && -1 != mesh_intervals() && 0 != mesh_intervals())
[ # # ][ # # ]
1107 : 0 : return measure() / mesh_intervals();
1108 : :
1109 : 0 : else return -1;
1110 : : }
1111 : :
1112 [ + - ][ + - ]: 156 : } // namespace MeshKit
1113 : :
|
