[MOAB-dev] commit/MOAB: iulian07: create the quad mesh for mbcslam, span between departure and arrival edges
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jul 9 15:18:31 CDT 2013
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/d018a2576a56/
Changeset: d018a2576a56
Branch: master
User: iulian07
Date: 2013-07-09 22:13:14
Summary: create the quad mesh for mbcslam, span between departure and arrival edges
it is not covering the initial domain, and it will be hard to make an
efficient intersection;
All quads will be more or less independent, and we cannot advance easily to the next quad, for advancing front intersection.
Affected #: 4 files
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 1166a7a..70d81ed 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -11,6 +11,7 @@
// right now, add a dependency to mbcoupler
#include "ElemUtil.hpp"
#include "moab/MergeMesh.hpp"
+#include "moab/ReadUtilIface.hpp"
// this is for sstream
#include <sstream>
@@ -991,5 +992,145 @@ ErrorCode enforce_convexity(Interface * mb, EntityHandle lset)
std::cout << brokenPolys << " concave polygons were decomposed in convex ones \n";
return MB_SUCCESS;
}
+ErrorCode create_span_quads(Interface * mb, EntityHandle euler_set, int rank)
+{
+ // first get all edges adjacent to polygons
+ Tag dpTag = 0;
+ std::string tag_name("DP");
+ ErrorCode rval = mb->tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE);
+ // if the tag does not exist, get out early
+ if (rval!=MB_SUCCESS)
+ return rval;
+ Range polygons;
+ rval = mb->get_entities_by_dimension(euler_set, 2, polygons);
+ if (MB_SUCCESS != rval)
+ return rval;
+ Range iniEdges;
+ rval = mb->get_adjacencies(polygons,
+ 1,
+ false,
+ iniEdges,
+ Interface::UNION);
+ if (MB_SUCCESS != rval)
+ return rval;
+ // now create some if missing
+ Range allEdges;
+ rval = mb->get_adjacencies(polygons,
+ 1,
+ true,
+ allEdges,
+ Interface::UNION);
+ // create the vertices at the DP points, and the quads after that
+ Range verts;
+ rval = mb->get_connectivity(polygons, verts);
+ if (MB_SUCCESS != rval)
+ return rval;
+ int num_verts = (int) verts.size();
+ // now see the departure points; to what boxes should we send them?
+ std::vector<double> dep_points(3*num_verts);
+ rval = mb->tag_get_data(dpTag, verts, (void*)&dep_points[0]);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // create vertices corresponding to dp locations
+ ReadUtilIface *read_iface;
+ rval = mb->query_interface(read_iface);
+ if (MB_SUCCESS != rval)
+ return rval;
+ std::vector<double *> coords;
+ EntityHandle start_vert, start_elem, *connect;
+ // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop over quads later
+ rval = read_iface->get_node_coords (3, num_verts, 0, start_vert, coords);
+ if (MB_SUCCESS != rval)
+ return rval;
+ // fill it up
+ for (int i=0; i<num_verts; i++)
+ {
+ // block from interleaved
+ coords[0][i]=dep_points[3*i];
+ coords[1][i]=dep_points[3*i+1];
+ coords[2][i]=dep_points[3*i+2];
+ }
+ // create quads; one quad for each edge
+ rval = read_iface->get_element_connect(allEdges.size(), 4, MBQUAD, 0, start_elem, connect);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ const EntityHandle * edge_conn = NULL;
+ int quad_index=0;
+ EntityHandle firstVertHandle = verts[0];// assume vertices are contiguous...
+ for (Range::iterator eit=allEdges.begin(); eit!=allEdges.end(); eit++, quad_index++)
+ {
+ EntityHandle edge=*eit;
+ int num_nodes;
+ rval = mb->get_connectivity(edge, edge_conn, num_nodes);
+ if (MB_SUCCESS != rval)
+ return rval;
+ connect[quad_index*4] = edge_conn[0];
+ connect[quad_index*4+1] = edge_conn[1];
+ // maybe some indexing in range?
+ connect[quad_index*4+2] = start_vert + edge_conn[1]- firstVertHandle;
+ connect[quad_index*4+3] = start_vert + edge_conn[0]- firstVertHandle;
+ }
+
+
+ Range quads(start_elem, start_elem+allEdges.size()-1);
+ EntityHandle outSet;
+ rval = mb->create_meshset(MESHSET_SET, outSet);
+ if (MB_SUCCESS != rval)
+ return rval;
+ mb->add_entities(outSet, quads);
+
+ Tag colTag;
+ rval = mb->tag_get_handle("COLOR_ID", 1, MB_TYPE_INTEGER,
+ colTag, MB_TAG_DENSE|MB_TAG_CREAT);
+ if (MB_SUCCESS != rval)
+ return rval;
+ int j=1;
+ for (Range::iterator itq=quads.begin(); itq!=quads.end(); itq++, j++)
+ {
+ EntityHandle q=*itq;
+ rval = mb->tag_set_data(colTag, &q, 1, &j);
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+ std::stringstream outf;
+ outf << "SpanQuads" << rank << ".h5m";
+ rval = mb->write_file(outf.str().c_str(), 0, 0, &outSet, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+ EntityHandle outSet2;
+ rval = mb->create_meshset(MESHSET_SET, outSet2);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ Range quadEdges;
+ rval = mb->get_adjacencies(quads,
+ 1,
+ true,
+ quadEdges,
+ Interface::UNION);
+ mb->add_entities(outSet2, quadEdges);
+
+ std::stringstream outf2;
+ outf2 << "SpanEdges" << rank << ".h5m";
+ rval = mb->write_file(outf2.str().c_str(), 0, 0, &outSet2, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+// maybe some clean up
+ mb->delete_entities(&outSet, 1);
+ mb->delete_entities(&outSet2, 1);
+ mb->delete_entities(quads);
+ Range new_edges=subtract(allEdges, iniEdges);
+ mb->delete_entities(new_edges);
+ new_edges=subtract(quadEdges, iniEdges);
+ mb->delete_entities(new_edges);
+ Range new_verts(start_vert, start_vert+num_verts);
+ mb->delete_entities(new_verts);
+
+ return MB_SUCCESS;
+
+}
} //namespace moab
diff --git a/tools/mbcslam/CslamUtils.hpp b/tools/mbcslam/CslamUtils.hpp
index c32a264..241fe84 100644
--- a/tools/mbcslam/CslamUtils.hpp
+++ b/tools/mbcslam/CslamUtils.hpp
@@ -130,5 +130,9 @@ void departure_point_case1(CartVect & arrival_point, double t, double delta_t, C
// maybe radius is not needed;
//
ErrorCode enforce_convexity(Interface * mb, EntityHandle set);
+
+// looking at DP tag, create the spanning quads, write a file (with rank) and
+// then delete the new entities (vertices) and the set of quads
+ErrorCode create_span_quads(Interface * mb, EntityHandle euler_set, int rank);
}
#endif /* CSLAMUTILS_HPP_ */
diff --git a/tools/mbcslam/Makefile.am b/tools/mbcslam/Makefile.am
index 49e54f6..68a1812 100644
--- a/tools/mbcslam/Makefile.am
+++ b/tools/mbcslam/Makefile.am
@@ -70,5 +70,7 @@ MOSTLYCLEANFILES = intersect1.h5m \
lagrIni.h5m \
intx.vtk \
spectral.vtk \
- intx1.vtk
+ intx1.vtk \
+ SpanEdges0.h5m \
+ SpanQuads0.h5m
diff --git a/tools/mbcslam/intx_mpas.cpp b/tools/mbcslam/intx_mpas.cpp
index 4bc742a..827871b 100644
--- a/tools/mbcslam/intx_mpas.cpp
+++ b/tools/mbcslam/intx_mpas.cpp
@@ -112,9 +112,9 @@ ErrorCode manufacture_lagrange_mesh_on_sphere(Interface * mb,
ptr_DP[2]=newPos[2];
ptr_DP+=3; // increment to the next node
}
-
return rval;
}
+
int main(int argc, char **argv)
{
@@ -170,6 +170,11 @@ int main(int argc, char **argv)
rval = manufacture_lagrange_mesh_on_sphere(&mb, euler_set);
if (MB_SUCCESS != rval)
return 1;
+ // create a set with quads corresponding to each initial edge spanned with the displacement field
+
+ rval = create_span_quads(&mb, euler_set, rank);
+ if (MB_SUCCESS != rval)
+ return 1;
EntityHandle covering_lagr_set;
rval = mb.create_meshset(MESHSET_SET, covering_lagr_set);
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list