[MOAB-dev] r1500 - in MOAB/trunk: . refiner

pebay at mcs.anl.gov pebay at mcs.anl.gov
Fri Dec 28 15:30:15 CST 2007


Author: pebay
Date: 2007-12-28 15:30:15 -0600 (Fri, 28 Dec 2007)
New Revision: 1500

Modified:
   MOAB/trunk/refiner/MBEdgeSizeEvaluator.cpp
   MOAB/trunk/refiner/MBEdgeSizeSimpleImplicit.cpp
   MOAB/trunk/refiner/MBEntityRefiner.cpp
   MOAB/trunk/refiner/MBEntityRefiner.hpp
   MOAB/trunk/refiner/MBSimplexTemplateRefiner.cpp
   MOAB/trunk/refiner/MBSimplexTemplateRefiner.hpp
   MOAB/trunk/test_mesh_refiner.cpp
Log:
ENH: implement and test the 2-simplex refinement


Modified: MOAB/trunk/refiner/MBEdgeSizeEvaluator.cpp
===================================================================
--- MOAB/trunk/refiner/MBEdgeSizeEvaluator.cpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBEdgeSizeEvaluator.cpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -96,5 +96,11 @@
                                                      const double* cm, void* tm, 
                                                      const double* c1, const void* t1 ) const
 {
+  (void)c0;
+  (void)t0;
+  (void)cm;
+  (void)tm;
+  (void)c1;
+  (void)t1;
 }
 

Modified: MOAB/trunk/refiner/MBEdgeSizeSimpleImplicit.cpp
===================================================================
--- MOAB/trunk/refiner/MBEdgeSizeSimpleImplicit.cpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBEdgeSizeSimpleImplicit.cpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -24,6 +24,9 @@
   double* p1, void* t1,
   const double* p2, const void* t2 )
 {
+  (void)t0;
+  (void)t1;
+  (void)t2;
   double L2 = 0.;
   double delta;
   int i;

Modified: MOAB/trunk/refiner/MBEntityRefiner.cpp
===================================================================
--- MOAB/trunk/refiner/MBEntityRefiner.cpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBEntityRefiner.cpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -8,6 +8,7 @@
 {  
   this->mesh = parentMesh;
   this->edge_size_evaluator = 0;
+  this->output_functor = 0;
   // By default, allow at most one subdivision per edge
   this->minimum_number_of_subdivisions = 0;
   this->maximum_number_of_subdivisions = 1;
@@ -18,6 +19,8 @@
 {
   if ( this->edge_size_evaluator )
     delete this->edge_size_evaluator;
+  if ( this->output_functor )
+    delete this->output_functor;
 }
 
 /**\fn bool MBEntityRefiner::refine_entity( MBEntityHandle )
@@ -146,11 +149,11 @@
 void MBEntityRefiner::update_heap_size()
 {
   unsigned long n = this->get_heap_size_bound( this->maximum_number_of_subdivisions );
-  this->coord_heap.reserve( 6 * n );
+  this->coord_heap.resize( 6 * n );
   if ( this->edge_size_evaluator )
     {
     unsigned long m = this->edge_size_evaluator->get_vertex_tag_size();
-    this->tag_heap.reserve( m * n );
+    this->tag_heap.resize( m * n );
     }
 }
 

Modified: MOAB/trunk/refiner/MBEntityRefiner.hpp
===================================================================
--- MOAB/trunk/refiner/MBEntityRefiner.hpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBEntityRefiner.hpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -69,6 +69,7 @@
 class MB_DLL_EXPORT MBEntityRefinerOutputFunctor
 {
 public:
+  virtual ~MBEntityRefinerOutputFunctor() { }
   virtual void operator () ( const double* vcoords, const void* vtags ) = 0;
   virtual void operator () ( MBEntityType etyp ) = 0;
 };

Modified: MOAB/trunk/refiner/MBSimplexTemplateRefiner.cpp
===================================================================
--- MOAB/trunk/refiner/MBSimplexTemplateRefiner.cpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBSimplexTemplateRefiner.cpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -6,6 +6,12 @@
 #include <iostream>
 #include <stack>
 
+// Static arrays holding parametric coordinates of element vertices
+static double MBVertexParametric[] = { 0., 0., 0. };
+static double MBEdgeParametric[]   = { 0., 0., 0.,   1., 0., 0. };
+static double MBTriParametric[]    = { 0., 0., 0.,   1., 0., 0.,   0., 1., 0. };
+static double MBTetParametric[]    = { 0., 0., 0.,   1., 0., 0.,   0., 1., 0.,   0., 0., 1. };
+
 #ifdef MB_DEBUG_TESSELLATOR
 #  define MB_TESSELLATOR_INCR_CASE_COUNT(cs) this->case_counts[cs]++
 #  define MB_TESSELLATOR_INCR_SUBCASE_COUNT(cs,sc) this->subcase_counts[cs][sc]++
@@ -36,28 +42,42 @@
     return false;
     }
   std::vector<double> entity_coords;
+  std::vector<char> entity_tags;
+  int ts = this->edge_size_evaluator->get_vertex_tag_size();
   entity_coords.resize( 6 * num_nodes );
+  if ( ts )
+    {
+    entity_tags.resize( num_nodes * ts );
+    }
+  MBEntityType etyp = this->mesh->type_from_handle( entity );
   // Have to make num_nodes calls to get_coords() because we need xyz interleaved with rst coords.
   for ( int n = 0; n < num_nodes; ++ n )
     {
-    if ( this->mesh->get_coords( &conn[n], 1, &entity_coords[3 * n + 3] ) != MB_SUCCESS )
+    if ( this->mesh->get_coords( &conn[n], 1, &entity_coords[6 * n + 3] ) != MB_SUCCESS )
       {
       return false;
       }
+    // Still need to get tags.
     }
-  // Still need to get tags.
 
-  switch ( this->mesh->type_from_handle( entity ) )
+  switch ( etyp )
     {
     case MBVERTEX:
-      this->refine_0_simplex( &entity_coords[0], 0 ); // FIXME
+      this->assign_parametric_coordinates( 1, MBVertexParametric, &entity_coords[0] );
+      this->refine_0_simplex( &entity_coords[0], &entity_tags[0] );
       rval = false;
       break;
     case MBEDGE:
-      rval = this->refine_1_simplex( 0,  &entity_coords[0], 0,  &entity_coords[6], 0 ); // FIXME
+      this->assign_parametric_coordinates( 2, MBEdgeParametric, &entity_coords[0] );
+      rval = this->refine_1_simplex( this->maximum_number_of_subdivisions,
+        &entity_coords[0], &entity_tags[0],  &entity_coords[6], &entity_tags[ts] );
       break;
     case MBTRI:
-      rval = this->refine_2_simplex( 0, 0,   0, 0,   0, 0,   0, 0 ); // FIXME
+      this->assign_parametric_coordinates( 3, MBTriParametric, &entity_coords[0] );
+      rval = this->refine_2_simplex( this->maximum_number_of_subdivisions, 7,
+        &entity_coords[ 0], &entity_tags[     0],
+        &entity_coords[ 6], &entity_tags[    ts],
+        &entity_coords[12], &entity_tags[2 * ts] );
       break;
     case MBQUAD:
       std::cerr << "Quadrilaterals not handled yet\n";
@@ -68,6 +88,7 @@
       rval = false;
       break;
     case MBTET:
+      this->assign_parametric_coordinates( 4, MBTetParametric, &entity_coords[0] );
       rval = this->refine_3_simplex( 0, 0, 0, 0, 0, 0, 0, 0, 0 ); // FIXME
       break;
     case MBPYRAMID:
@@ -184,6 +205,12 @@
   if ( max_depth-- > 0 )
     {
     int i;
+    midpt0c = this->heap_coord_storage();
+    midpt1c = this->heap_coord_storage();
+    midpt2c = this->heap_coord_storage();
+    midpt0t = this->heap_tag_storage();
+    midpt1t = this->heap_tag_storage();
+    midpt2t = this->heap_tag_storage();
     for ( i = 0; i < 6; ++i )
       {
       midpt0c[i] = ( v0[i] + v1[i] ) / 2.;
@@ -1419,6 +1446,15 @@
   return true;
 }
 
+/**\brief This is used by refine_entity to assign parametric coordinates to corners of each element.
+  */
+void MBSimplexTemplateRefiner::assign_parametric_coordinates( int num_nodes, const double* src, double* tgt )
+{
+  for ( int i = 0; i < num_nodes; ++i, src +=3, tgt += 6 )
+    for ( int j = 0; j < 3; ++j )
+      tgt[j] = src[j];
+}
+
 /**\brief Returns true if || a0a1 || < || b0b1 ||
   *
   * We use this to test which triangulation has the best

Modified: MOAB/trunk/refiner/MBSimplexTemplateRefiner.hpp
===================================================================
--- MOAB/trunk/refiner/MBSimplexTemplateRefiner.hpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/refiner/MBSimplexTemplateRefiner.hpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -60,13 +60,14 @@
                          double* v1, void* t1, 
                          double* v2, void* t2,
                          double* v3, void* t3 );
-  static bool compare_Hopf_cross_string_dist( const double* v00, const double* v01, const double* v10, const double* v11 );
   void evaluate_tags_at_facepoint( const double* c0, const void* t0,
 				   const double* c1, const void* t1,
 				   const double* c2, const void* t2,
 				   const double* cm, void* tm ) const;
 
   int best_tets( int* alternates, double*[14], int, int ) { return alternates[0]; }
+  void assign_parametric_coordinates( int num_nodes, const double* src, double* tgt );
+  static bool compare_Hopf_cross_string_dist( const double* v00, const double* v01, const double* v10, const double* v11 );
 };
 #endif // MB_SIMPLEXTEMPLATEREFINER_H
 

Modified: MOAB/trunk/test_mesh_refiner.cpp
===================================================================
--- MOAB/trunk/test_mesh_refiner.cpp	2007-12-28 20:58:04 UTC (rev 1499)
+++ MOAB/trunk/test_mesh_refiner.cpp	2007-12-28 21:30:15 UTC (rev 1500)
@@ -1,25 +1,62 @@
 #include "MBCore.hpp"
 #include "MBEdgeSizeSimpleImplicit.hpp"
+#include "MBSimplexTemplateRefiner.hpp"
 #include "MBInterface.hpp"
 
+#include <iostream>
+
+class MBTestOutputFunctor : public MBEntityRefinerOutputFunctor
+{
+  virtual void operator () ( const double* vcoords, const void* vtags )
+    {
+    std::cout << "[ " << vcoords[0];
+    for ( int i = 1; i < 6; ++i )
+      std::cout << ", " << vcoords[i];
+    std::cout << " ]\n";
+    }
+
+  virtual void operator () ( MBEntityType etyp )
+    {
+    std::cout << "---------- " << etyp << "\n\n";
+    }
+};
+
 int TestMeshRefiner( int argc, char* argv[] )
 {
   MBInterface* iface = new MBCore;
-  MBEdgeSizeSimpleImplicit eval( iface );
-  double p0[6] = { 0.1, 0.0, 0.0,  0.1, 0.0, 0.0 };
-  double p1[6] = { 0.6, 0.0, 0.0,  0.6, 0.0, 0.0 };
-  double p2[6] = { 1.1, 0.0, 0.0,  1.1, 0.0, 0.0 };
-  if ( eval.evaluate_edge( p0, 0, p1, 0, p2, 0 ) )
+  MBEdgeSizeSimpleImplicit* eval = new MBEdgeSizeSimpleImplicit( iface );
+
+  double p0[6] = { 0.0, 0.0, 0.0,  0.0, 0.0, 0.0 };
+  double p1[6] = { 0.5, 0.0, 0.0,  0.5, 0.0, 0.0 };
+  double p2[6] = { 1.0, 0.0, 0.0,  1.0, 0.0, 0.0 };
+  double p3[6] = { 0.6, 2.0, 0.0,  0.6, 2.0, 0.0 };
+
+  if ( eval->evaluate_edge( p0, 0, p1, 0, p2, 0 ) )
     {
     return 1;
     }
 
-  eval.set_ratio( 2. );
-  if ( ! eval.evaluate_edge( p0, 0, p1, 0, p2, 0 ) )
+  eval->set_ratio( 2. );
+  if ( ! eval->evaluate_edge( p0, 0, p1, 0, p2, 0 ) )
     {
     return 1;
     }
 
+  MBEntityHandle node_handles[3];
+  MBEntityHandle tri_handle;
+
+  iface->create_vertex( p0, node_handles[0] );
+  iface->create_vertex( p2, node_handles[1] );
+  iface->create_vertex( p3, node_handles[2] );
+
+  iface->create_element( MBTRI, node_handles, 3, tri_handle );
+
+  MBSimplexTemplateRefiner eref( iface );
+  MBTestOutputFunctor* ofunc = new MBTestOutputFunctor;
+  eref.set_edge_size_evaluator( eval );
+  eref.set_output_functor( ofunc );
+  eref.refine_entity( tri_handle );
+
   delete iface;
   return 0;
 }




More information about the moab-dev mailing list