[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