Merge pull request #4074 from roystgnr/polygon_elements

(Arbitrary-) Polygon elements

Author: roystgnr

Makefile.in

@@ -744,9 +744,11 @@ include_HEADERS = \
         geom/elem_range.h \
         geom/elem_side_builder.h \
         geom/face.h \
+        geom/face_c0polygon.h \
         geom/face_inf_quad.h \
         geom/face_inf_quad4.h \
         geom/face_inf_quad6.h \
+        geom/face_polygon.h \
         geom/face_quad.h \
         geom/face_quad4.h \
         geom/face_quad4_shell.h \

include/Makefile.in

@@ -79,6 +79,8 @@ enum ElemType : int {
                PYRAMID18 = 37,
                // Another shell elem
                QUADSHELL9 = 38,
+               // An arbitrary polygon with N EDGE2 sides
+               C0POLYGON = 39,
                // Invalid
                INVALID_ELEM};   // should always be last
 

include/enums/enum_elem_type.h

@@ -488,12 +488,25 @@ class FEAbstract : public ReferenceCountedObject<FEAbstract>
    */
   virtual unsigned int n_quadrature_points () const = 0;
 
+  /**
+   * \returns The element that the current shape functions have been
+   * calculated for.  Useful in determining when shape functions must
+   * be recomputed.
+   */
+  const Elem * get_elem()  const { return _elem; }
+
   /**
    * \returns The element type that the current shape functions
-   * have been calculated for.  Useful in determining when shape
-   * functions must be recomputed.
+   * have been calculated for, or \p INVALID_ELEM if no such element
+   * exists.  Useful in determining when shape functions must be
+   * recomputed.
+   *
+   * This is generally redundant with _elem->type(), but must be
+   * cached separately for cases (such as internal FE use in
+   * QComposite) where _elem might be a dangling pointer to a
+   * temporary.
    */
-  ElemType get_type()  const { return elem_type; }
+  ElemType get_type() const { return _elem_type; }
 
   /**
    * \returns The p refinement level that the current shape
@@ -717,10 +730,14 @@ class FEAbstract : public ReferenceCountedObject<FEAbstract>
   FEType fe_type;
 
   /**
-   * The element type the current data structures are
-   * set up for.
+   * The element type the current data structures were set up for.
+   */
+  ElemType _elem_type;
+
+  /**
+   * The element the current data structures were set up for.
    */
-  ElemType elem_type;
+  const Elem * _elem;
 
   /**
    * The element p-refinement level the current data structures are

include/fe/fe_abstract.h

@@ -237,13 +237,25 @@ class Elem : public ReferenceCountedObject<Elem>,
    */
   static const subdomain_id_type invalid_subdomain_id;
 
+  /**
+   * \returns true iff this element type can vary in topology (e.g.
+   * have different numbers of sides and/or nodes) at runtime.  For
+   * such general polygons or polyhedra, APIs which assume a fixed
+   * topology are not safe to use.
+   */
+  virtual bool runtime_topology() const { return false; }
+
   /**
    * \returns A pointer to the "reference element" associated
    * with this element.  The reference element is the image of this
    * element in reference parametric space. Importantly, it is *not*
    * an actual element in the mesh, but rather a Singleton-type
    * object, so for example all \p Quad4 elements share the same
    * \p reference_elem().
+   *
+   * If the element is of a type that can admit multiple topologies,
+   * such as a Polygon subtype, then there is no reference element;
+   * for such types this method should not be used.
    */
   const Elem * reference_elem () const;
 
@@ -608,6 +620,10 @@ class Elem : public ReferenceCountedObject<Elem>,
   /**
    * This array maps the integer representation of the \p ElemType enum
    * to the number of nodes in the element.
+   *
+   * This is only usable for simple types for which the node number
+   * is fixed; for more general types like Polygon subclasses an actual
+   * instantiated Elem must be queried.
    */
   static const unsigned int type_to_n_nodes_map[INVALID_ELEM];
 
@@ -639,6 +655,10 @@ class Elem : public ReferenceCountedObject<Elem>,
   /**
    * This array maps the integer representation of the \p ElemType enum
    * to the number of sides on the element.
+   *
+   * This is only usable for simple types for which the node number
+   * is fixed; for more general types like Polygon subclasses an actual
+   * instantiated Elem must be queried.
    */
   static const unsigned int type_to_n_sides_map[INVALID_ELEM];
 
@@ -693,6 +713,10 @@ class Elem : public ReferenceCountedObject<Elem>,
   /**
    * This array maps the integer representation of the \p ElemType enum
    * to the number of edges on the element.
+   *
+   * This is only usable for simple types for which the node number
+   * is fixed; for more general types like Polygon subclasses an actual
+   * instantiated Elem must be queried.
    */
   static const unsigned int type_to_n_edges_map[INVALID_ELEM];
 
@@ -1057,8 +1081,8 @@ class Elem : public ReferenceCountedObject<Elem>,
   { libmesh_not_implemented(); return std::make_pair(0.,0.); }
 
   /**
-   * \returns \p true if the point p is contained in this element,
-   * false otherwise.
+   * \returns \p true if the physical point p is contained in this
+   * element, false otherwise.
    *
    * For linear elements, performs an initial tight bounding box check
    * (as an optimization step) and (if that passes) then uses the
@@ -1844,14 +1868,12 @@ class Elem : public ReferenceCountedObject<Elem>,
 
 #endif
 
-
-
-
   /**
    * \returns An Elem of type \p type wrapped in a smart pointer.
    */
   static std::unique_ptr<Elem> build (const ElemType type,
                                       Elem * p=nullptr);
+
   /**
    * Calls the build() method above with a nullptr parent, and
    * additionally sets the newly-created Elem's id. This can be useful
@@ -1860,6 +1882,20 @@ class Elem : public ReferenceCountedObject<Elem>,
   static std::unique_ptr<Elem> build_with_id (const ElemType type,
                                               dof_id_type id);
 
+  /**
+   * \returns An Elem of the same type as \p this, wrapped in a smart
+   * pointer.
+   *
+   * This is not a complete clone() method (since e.g. it does not set
+   * node pointers; the standard use case reassigns node pointers from
+   * a different mesh), but it is necessary to use this instead of
+   * build() for runtime-polymorphic elements like Polygon subtypes
+   * whose "type" depends on more than their type(), and it is useful
+   * to use this for elements whose id, unique_id, extra integers,
+   * etc. should be preserved in the near-clone.
+   */
+  virtual std::unique_ptr<Elem> disconnected_clone () const;
+
   /**
    * Returns the number of independent permutations of element nodes -
    * e.g. a cube can be reoriented to put side 0 where side N is (for
@@ -2310,7 +2346,8 @@ Elem::Elem(const unsigned int nn,
   // If this ever legitimately fails we need to increase max_n_nodes
   libmesh_assert_less_equal(nn, max_n_nodes);
 
-  // Initialize the nodes data structure
+  // Initialize the nodes data structure if we're given a pointer to
+  // memory for it.
   if (_nodes)
     {
       for (unsigned int n=0; n<nn; n++)
@@ -2320,27 +2357,31 @@ Elem::Elem(const unsigned int nn,
   // Initialize the neighbors/parent data structure
   // _elemlinks = new Elem *[ns+1];
 
-  // We now require that we get allocated data from a subclass
-  libmesh_assert (_elemlinks);
-
-  _elemlinks[0] = p;
+  // Initialize the elements data structure if we're given a pointer
+  // to memory for it.  If we *weren't* given memory for it, e.g.
+  // because a subclass like an arbitrary Polygon needs to
+  // heap-allocate this memory, then that subclass will have to handle
+  // this initialization too.
+  if (_elemlinks)
+    {
+      _elemlinks[0] = p;
 
-  for (unsigned int n=1; n<ns+1; n++)
-    _elemlinks[n] = nullptr;
+      for (unsigned int n=1; n<ns+1; n++)
+        _elemlinks[n] = nullptr;
 
-  // Optionally initialize data from the parent
-  if (this->parent() != nullptr)
-    {
-      this->subdomain_id() = this->parent()->subdomain_id();
-      this->processor_id() = this->parent()->processor_id();
-      _map_type = this->parent()->_map_type;
-      _map_data = this->parent()->_map_data;
-    }
+      // Optionally initialize data from the parent
+      if (this->parent())
+        {
+          this->subdomain_id() = this->parent()->subdomain_id();
+          this->processor_id() = this->parent()->processor_id();
+          _map_type = this->parent()->_map_type;
+          _map_data = this->parent()->_map_data;
 
 #ifdef LIBMESH_ENABLE_AMR
-  if (this->parent())
-    this->set_p_level(this->parent()->p_level());
+          this->set_p_level(this->parent()->p_level());
 #endif
+        }
+    }
 }