[MOAB-dev] r1418 - MOAB/trunk
tautges at mcs.anl.gov
tautges at mcs.anl.gov
Wed Nov 21 09:22:26 CST 2007
Author: tautges
Date: 2007-11-21 09:22:26 -0600 (Wed, 21 Nov 2007)
New Revision: 1418
Modified:
MOAB/trunk/DualTool.cpp
MOAB/trunk/DualTool.hpp
MOAB/trunk/MeshTopoUtil.cpp
MOAB/trunk/MeshTopoUtil.hpp
Log:
MeshTopoUtil: in split_entities_manifold, provide option to tell which higher-dimensional entity the new entity goes with
DualTool: forgot to include faces and hexes in entities to update when splitting node nonmanifold
Modified: MOAB/trunk/DualTool.cpp
===================================================================
--- MOAB/trunk/DualTool.cpp 2007-11-21 13:17:58 UTC (rev 1417)
+++ MOAB/trunk/DualTool.cpp 2007-11-21 15:22:26 UTC (rev 1418)
@@ -1693,8 +1693,9 @@
// reconstruct dual
result = construct_hex_dual(hexes);
if (MB_SUCCESS != result) return result;
+
+ return check_dual_adjs();
- return MB_SUCCESS;
}
MBErrorCode DualTool::foc_get_ents(MBEntityHandle ocl,
@@ -1830,8 +1831,22 @@
//=============== split faces
+ for (int i = 0; i < 2; i++) {
+ // get a hex in star_dp2[0] that's adj to this split quad, to tell
+ // mtu which one the orig should go with
+ MBEntityHandle gowith_hex = 0;
+ for (std::vector<MBEntityHandle>::iterator vit = star_dp2[0].begin();
+ vit != star_dp2[0].end(); vit++) {
+ if (mtu.common_entity(*vit, split_quads[i], 2)) {
+ gowith_hex = *vit;
+ break;
+ }
+ }
+
// split manifold each of the split_quads, and put the results on the merge list
- result = mtu.split_entities_manifold(split_quads, 2, new_quads, NULL); RR;
+ result = mtu.split_entities_manifold(split_quads+i, 1, new_quads+i, NULL,
+ (gowith_hex ? &gowith_hex : NULL)); RR;
+ }
// make ranges of faces which need to be explicitly adj to old, new
// edge; faces come from stars and new_quads (which weren't in the stars)
@@ -1873,8 +1888,6 @@
//=============== split node
- //=============== prepare for splitting 2 edges/1 node part
-
// if we're splitting 2 edges, there might be other edges that have the split
// node; also need to know which side they're on
MBRange addl_edges[2];
@@ -1888,6 +1901,15 @@
addl_edges[1].insert(split_edges[i]);
}
+ // same for star faces and hexes
+ for (int i = 0; i < 2; i++) {
+ std::copy(star_dp1[i].begin(), star_dp1[i].end(), mb_range_inserter(addl_edges[i]));
+ std::copy(star_dp2[i].begin(), star_dp2[i].end(), mb_range_inserter(addl_edges[i]));
+ }
+
+ // finally, new quads
+ for (int i = 0; i < 2; i++) addl_edges[0].insert(new_quads[i]);
+
// now split the node too
result = mtu.split_entity_nonmanifold(split_node, addl_edges[1],
addl_edges[0], new_node); RR;
@@ -2877,3 +2899,67 @@
return MB_SUCCESS;
}
+
+MBErrorCode DualTool::check_dual_adjs()
+{
+ // check primal-dual correspondence
+
+ // get the primal entities
+ MBRange pents[4];
+ MBErrorCode result = mbImpl->get_entities_by_type(0, MBHEX, pents[3]); RR;
+ for (int i = 2; i >= 0; i--) {
+ result = mbImpl->get_adjacencies(pents[3], 2, false, pents[2],
+ MBInterface::UNION); RR;
+ }
+
+ // for each primal entity of dimension pd
+#define PRENT(ent) MBCN::EntityTypeName(TYPE_FROM_HANDLE(ent)) << " " \
+ << ID_FROM_HANDLE(ent)
+ MBErrorCode overall_result = MB_SUCCESS;
+ for (int pd = 1; pd <= 3; pd++) {
+ for (MBRange::iterator prit = pents[pd].begin(); prit != pents[pd].end(); prit++) {
+ // get corresponding dual entity of dimension dd = 3-pd
+ MBEntityHandle dual_ent = get_dual_entity(*prit);
+ if (0 == dual_ent)
+ std::cerr << "Problem getting dual entity for " << PRENT(*prit) << std::endl;
+
+ // for each sub dimension sd = 0..pd-1
+ for (int sd = 0; sd < pd; sd++) {
+ MBRange R1, R2, R3;
+ // R1 = entities bounding primal entity of dim sd
+ result = mbImpl->get_adjacencies(&(*prit), 1, sd, false, R1); RR;
+
+ // R2 = entities bounded by dual entity, of dim 3-sd
+ result = mbImpl->get_adjacencies(&dual_ent, 1, 3-sd, false, R2); RR;
+
+
+ if (R1.size() != R2.size()) {
+ std::cerr << PRENT(*prit) << ": number of adj ents in "
+ << "primal/dual don't agree for dimension " << sd << "." << std::endl;
+ overall_result = MB_FAILURE;
+ }
+
+ // for each entity in R1, get its dual and look for it in R2
+ for (MBRange::iterator r1it = R1.begin(); r1it != R1.end(); r1it++) {
+ MBEntityHandle tmp_dual = get_dual_entity(*r1it);
+ if (R2.find(tmp_dual) == R2.end()) {
+ std::cerr << PRENT(*prit) << ": adj entity " << PRENT(*r1it)
+ << " isn't adjacent in dual." << std::endl;
+ overall_result = MB_FAILURE;
+ }
+ }
+ // ditto for R2
+ for (MBRange::iterator r2it = R2.begin(); r2it != R2.end(); r2it++) {
+ MBEntityHandle tmp_prim = get_dual_entity(*r2it);
+ if (R1.find(tmp_prim) == R1.end()) {
+ std::cerr << PRENT(*prit) << ": adj entity " << PRENT(*r2it)
+ << " isn't adjacent in primal." << std::endl;
+ overall_result = MB_FAILURE;
+ }
+ }
+ }
+ }
+ }
+
+ return overall_result;
+}
Modified: MOAB/trunk/DualTool.hpp
===================================================================
--- MOAB/trunk/DualTool.hpp 2007-11-21 13:17:58 UTC (rev 1417)
+++ MOAB/trunk/DualTool.hpp 2007-11-21 15:22:26 UTC (rev 1418)
@@ -217,6 +217,9 @@
//! delete all the dual data
MBErrorCode delete_whole_dual();
+ //! check dual-primal adjacencies
+ MBErrorCode check_dual_adjs();
+
private:
//! construct dual vertices for specified regions
Modified: MOAB/trunk/MeshTopoUtil.cpp
===================================================================
--- MOAB/trunk/MeshTopoUtil.cpp 2007-11-21 13:17:58 UTC (rev 1417)
+++ MOAB/trunk/MeshTopoUtil.cpp 2007-11-21 15:22:26 UTC (rev 1418)
@@ -535,7 +535,8 @@
MBErrorCode MeshTopoUtil::split_entities_manifold(MBEntityHandle *entities,
const int num_entities,
MBEntityHandle *new_entities,
- MBRange *fill_entities)
+ MBRange *fill_entities,
+ MBEntityHandle *gowith_ents)
{
// split entities by duplicating them; splitting manifold means that there is at
// most two higher-dimension entities bounded by a given entity; after split, the
@@ -591,17 +592,27 @@
}
else {
- // first the new entity; if there's only one up_adj of this dimension, make the
- // new entity adjacent to it, so that any bdy information is preserved on the
- // original entity
- tmp_result = mbImpl->remove_adjacencies(entities[i], &(*(up_adjs[dim].begin())), 1);
+ // get the two up-elements
+ MBEntityHandle up_elem1 = *(up_adjs[dim].begin()),
+ up_elem2 = (up_adjs[dim].size() > 1 ? *(up_adjs[dim].rbegin()) : 0);
+
+ // if two, and a gowith entity was input, make sure the new entity goes with
+ // that one
+ if (gowith_ents && up_elem2 &&
+ gowith_ents[i] != up_elem1 && gowith_ents[i] == up_elem2) {
+ MBEntityHandle tmp_elem = up_elem1;
+ up_elem1 = up_elem2;
+ up_elem2 = tmp_elem;
+ }
+
+ tmp_result = mbImpl->remove_adjacencies(entities[i], &up_elem1, 1);
// (ok if there's an error, that just means there wasn't an explicit adj)
- MBEntityHandle tmp_ent = *(up_adjs[dim].begin());
- tmp_result = mbImpl->add_adjacencies(new_entity, &tmp_ent, 1, false); TC;
- if (up_adjs[dim].size() < 2) continue;
+
+ tmp_result = mbImpl->add_adjacencies(new_entity, &up_elem1, 1, false); TC;
+ if (!up_elem2) continue;
+
// add adj to other up_adj
- tmp_ent = *(up_adjs[dim].rbegin());
- tmp_result = mbImpl->add_adjacencies(entities[i], &tmp_ent, 1, false); TC;
+ tmp_result = mbImpl->add_adjacencies(entities[i], &up_elem2, 1, false); TC;
}
}
Modified: MOAB/trunk/MeshTopoUtil.hpp
===================================================================
--- MOAB/trunk/MeshTopoUtil.hpp 2007-11-21 13:17:58 UTC (rev 1417)
+++ MOAB/trunk/MeshTopoUtil.hpp 2007-11-21 15:22:26 UTC (rev 1418)
@@ -138,11 +138,17 @@
\param new_entities New entities, in order of correspondence to that of entities
\param fill_entities If non-NULL, create an entity of next higher dimension to fill the gap,
passing it back in *fill_entities
+ \param gowith_ents If non-NULL, each of the new entities will adj to the
+ corresponding gowith entities after the split; this parameter is
+ ignored for boundary split entities; in that case, the split entity
+ remains on the boundary (i.e. not adj to any entity of higher
+ dimension). Dimension of gowith_ents must be the same as entities.
*/
MBErrorCode split_entities_manifold(MBEntityHandle *entities,
const int num_entities,
MBEntityHandle *new_entities,
- MBRange *fill_entities);
+ MBRange *fill_entities,
+ MBEntityHandle *gowith_ents = NULL);
//! return whether entity is equivalent to any other of same type and same vertices;
//! if equivalent entity is found, it's returned in equiv_ents and return value is true,
More information about the moab-dev
mailing list