elem.c

一个用来实现偏微分方程中网格的计算库

          for (unsigned int s=0; s != child->n_sides(); ++s)
            {
	      Elem *child_neigh = child->neighbor(s);
              if (child_neigh &&
		  (child_neigh == neighbor ||
		   (child_neigh->parent() == neighbor &&
		    child_neigh->is_ancestor_of(subneighbor))))
                child->family_tree_by_subneighbor (family, child_neigh,
                                                   subneighbor, false);
            }
      }
}



void Elem::active_family_tree_by_neighbor (std::vector<const Elem*>& family,
                                           const Elem* neighbor,
			                   const bool reset) const
{
  // Clear the vector if the flag reset tells us to.
  if (reset)
    family.clear();

  // This only makes sense if we're already a neighbor
  if (this->level() >= neighbor->level())
    libmesh_assert (this->has_neighbor(neighbor));

  // Add an active element to the family tree.
  if (this->active())
    family.push_back(this);

  // Or recurse into an ancestor element's children.
  // Do not clear the vector any more.
  else if (this->has_children())
    for (unsigned int c=0; c<this->n_children(); c++)
      {
        Elem *child = this->child(c);
        if (child != remote_elem && child->has_neighbor(neighbor))
          child->active_family_tree_by_neighbor (family, neighbor, false);
      }
}



unsigned int Elem::min_p_level_by_neighbor(const Elem* neighbor,
                                           unsigned int current_min) const
{
  libmesh_assert(!this->subactive());
  libmesh_assert(neighbor->active());

  // If we're an active element this is simple
  if (this->active())
    return std::min(current_min, this->p_level());

  libmesh_assert(has_neighbor(neighbor));

  // The p_level() of an ancestor element is already the minimum
  // p_level() of its children - so if that's high enough, we don't
  // need to examine any children.
  if (current_min <= this->p_level())
    return current_min;

  unsigned int min_p_level = current_min;

  for (unsigned int c=0; c<this->n_children(); c++)
    {
      const Elem* const child = this->child(c);
      if (child != remote_elem && child->has_neighbor(neighbor))
        min_p_level =
	  child->min_p_level_by_neighbor(neighbor,
                                         min_p_level);
    }

  return min_p_level;
}


unsigned int Elem::min_new_p_level_by_neighbor(const Elem* neighbor,
                                               unsigned int current_min) const
{
  libmesh_assert(!this->subactive());
  libmesh_assert(neighbor->active());

  // If we're an active element this is simple
  if (this->active())
    {
      unsigned int new_p_level = this->p_level();
      if (this->p_refinement_flag() == Elem::REFINE)
        new_p_level += 1;
      if (this->p_refinement_flag() == Elem::COARSEN)
        {
          libmesh_assert (new_p_level > 0);
          new_p_level -= 1;
        }
      return std::min(current_min, new_p_level);
    }

  libmesh_assert(has_neighbor(neighbor));

  unsigned int min_p_level = current_min;

  for (unsigned int c=0; c<this->n_children(); c++)
    {
      const Elem* const child = this->child(c);
      if (child && child != remote_elem)
        if (child->has_neighbor(neighbor))
          min_p_level =
	    child->min_new_p_level_by_neighbor(neighbor,
                                               min_p_level);
    }

  return min_p_level;
}

#endif // #ifdef ENABLE_AMR



bool Elem::contains_point (const Point& p) const
{
  // Declare a basic FEType.  Will ue a Lagrange
  // element by default.
  FEType fe_type(this->default_order());
  
  const Point mapped_point = FEInterface::inverse_map(this->dim(),
						      fe_type,
						      this,
						      p,
						      1.e-4,
						      false);

  return FEInterface::on_reference_element(mapped_point, this->type());
}



void Elem::nullify_neighbors ()
{
  // Tell any of my neighbors about my death...
  // Looks strange, huh?
  for (unsigned int n=0; n<this->n_neighbors(); n++)
    {
      Elem* neighbor = this->neighbor(n);
      if (neighbor && neighbor != remote_elem)
        {
	  // Note:  it is possible that I see the neighbor
	  // (which is coarser than me)
	  // but they don't see me, so avoid that case.
	  if (neighbor->level() == this->level())
	    {	
	      const unsigned int w_n_a_i = neighbor->which_neighbor_am_i(this);
              libmesh_assert (w_n_a_i < neighbor->n_neighbors());
	      neighbor->set_neighbor(w_n_a_i, NULL);
	      this->set_neighbor(n, NULL);
	    }
        }
    }
}




unsigned int Elem::n_second_order_adjacent_vertices (const unsigned int) const
{
  // for linear elements, always return 0
  return 0;
}



unsigned short int Elem::second_order_adjacent_vertex (const unsigned int,
						       const unsigned int) const
{
  // for linear elements, always return 0
  return 0;
}



std::pair<unsigned short int, unsigned short int>
Elem::second_order_child_vertex (const unsigned int) const
{
  // for linear elements, always return 0
  return std::pair<unsigned short int, unsigned short int>(0,0);
}



ElemType Elem::first_order_equivalent_type (const ElemType et)
{ 
  switch (et)
    {
    case EDGE2:
    case EDGE3:
    case EDGE4:
      return EDGE2;
    case TRI3:
    case TRI6:
      return TRI3;
    case QUAD4:
    case QUAD8:
    case QUAD9:
      return QUAD4;
    case TET4:
    case TET10:
      return TET4;
    case HEX8:
    case HEX27:
    case HEX20:
      return HEX8;
    case PRISM6:
    case PRISM15:
    case PRISM18:
      return PRISM6;

#ifdef ENABLE_INFINITE_ELEMENTS

    case INFQUAD4:
    case INFQUAD6:
      return INFQUAD4;
    case INFHEX8:
    case INFHEX16:
    case INFHEX18:
      return INFHEX8;
    case INFPRISM6:
    case INFPRISM12:
      return INFPRISM6;

#endif

    default:
      // unknown element
      return INVALID_ELEM;
    }
}



ElemType Elem::second_order_equivalent_type (const ElemType et,
					     const bool full_ordered)
{ 
  /* for second-order elements, always return \p INVALID_ELEM
   * since second-order elements should not be converted
   * into something else.  Only linear elements should 
   * return something sensible here
   */
  switch (et)
    {
    case EDGE2:
      {
	// full_ordered not relevant
	return EDGE3;
      }

    case TRI3:
      {
	// full_ordered not relevant
	return TRI6;
      }

    case QUAD4:
      {
	if (full_ordered)
	  return QUAD9;
	else 
	  return QUAD8;
      }

    case TET4:
      {
	// full_ordered not relevant
	return TET10;
      }

    case HEX8:
      {
	// see below how this correlates with INFHEX8
	if (full_ordered)
	  return HEX27;
	else 
	  return HEX20;
      }

    case PRISM6:
      {
	if (full_ordered)
	  return PRISM18;
	else 
	  return PRISM15;
      }

    case PYRAMID5:
      {
	// libmesh_error(); 
	return INVALID_ELEM;
      }



#ifdef ENABLE_INFINITE_ELEMENTS

    // infinite elements
    case INFEDGE2:
      {
	// libmesh_error(); 
	return INVALID_ELEM;
      }

    case INFQUAD4:
      {
	// full_ordered not relevant
	return INFQUAD6;
      }

    case INFHEX8:
      {
	/*
	 * Note that this matches with \p Hex8:
	 * For full-ordered, \p InfHex18 and \p Hex27
	 * belong together, and for not full-ordered,
	 * \p InfHex16 and \p Hex20 belong together.
	 */
	if (full_ordered)
	  return INFHEX18;
	else 
	  return INFHEX16;
      }

    case INFPRISM6:
      {
	// full_ordered not relevant
	return INFPRISM12;
      }

#endif


    default:
      {
	// second-order element
	return INVALID_ELEM;
      }
    }
}



Elem::side_iterator Elem::boundary_sides_begin()
{
  Predicates::BoundarySide<SideIter> bsp;
  return side_iterator(this->_first_side(), this->_last_side(), bsp);
}




Elem::side_iterator Elem::boundary_sides_end()
{
  Predicates::BoundarySide<SideIter> bsp;
  return side_iterator(this->_last_side(), this->_last_side(), bsp);
}




Real Elem::volume () const
{
  // The default implementation builds a finite element of the correct
  // order and sums up the JxW contributions.  This can be expensive,
  // so the various element types can overload this method and compute
  // the volume more efficiently.
  FEType fe_type (this->default_order() , LAGRANGE);

  AutoPtr<FEBase> fe (FEBase::build(this->dim(),
				    fe_type));

   const std::vector<Real>& JxW = fe->get_JxW();
   
  // The default quadrature rule should integrate the mass matrix,
  // thus it should be plenty to compute the area
  QGauss qrule (this->dim(), fe_type.default_quadrature_order());

  fe->attach_quadrature_rule(&qrule);

  fe->reinit(this);

  Real vol=0.;
  for (unsigned int qp=0; qp<qrule.n_points(); ++qp)
    vol += JxW[qp];
  
  return vol;
  
}


// ------------------------------------------------------------
// Elem::PackedElem static data
const unsigned int Elem::PackedElem::header_size = 10;


// Elem::PackedElem member funcions
void Elem::PackedElem::pack (std::vector<int> &conn, const Elem* elem)
{
  libmesh_assert (elem != NULL);
  
  // we can do at least this good. note that hopefully in general
  // the user will already have reserved the full space, which will render
  // this redundant
  conn.reserve (conn.size() + Elem::PackedElem::header_size + elem->n_nodes());

#ifdef ENABLE_AMR
  conn.push_back (static_cast<int>(elem->level()));
  conn.push_back (static_cast<int>(elem->p_level()));
  conn.push_back (static_cast<int>(elem->refinement_flag()));
  conn.push_back (static_cast<int>(elem->p_refinement_flag()));
#else
  conn.push_back (0);
  conn.push_back (0);
  conn.push_back (0);
  conn.push_back (0);
#endif
  conn.push_back (static_cast<int>(elem->type()));
  conn.push_back (static_cast<int>(elem->processor_id()));
  conn.push_back (static_cast<int>(elem->subdomain_id()));
  conn.push_back (elem->id());
		
#ifdef ENABLE_AMR
  // use parent_ID of -1 to indicate a level 0 element
  if (elem->level() == 0)
    {
      conn.push_back(-1);
      conn.push_back(-1);
    }
  else
    {
      conn.push_back(elem->parent()->id());
      conn.push_back(elem->parent()->which_child_am_i(elem));
    }
#else
  conn.push_back (-1);
  conn.push_back (-1);
#endif
  
  for (unsigned int n=0; n<elem->n_nodes(); n++)
    conn.push_back (elem->node(n));		
}



Elem * Elem::PackedElem::unpack (MeshBase &mesh, Elem *parent) const
{  
  
  Elem *elem = Elem::build(this->type(),parent).release();
  libmesh_assert (elem);

#ifdef ENABLE_AMR
  if (this->level() != 0) 
    {
      libmesh_assert (parent != NULL);
      parent->add_child(elem, this->which_child_am_i());
      libmesh_assert (parent->type() == elem->type());
      libmesh_assert (parent->child(this->which_child_am_i()) == elem);
    }
#endif

  // Assign the IDs
#ifdef ENABLE_AMR
  elem->set_p_level(this->p_level());
  elem->set_refinement_flag(this->refinement_flag());
  elem->set_p_refinement_flag(this->p_refinement_flag());
  libmesh_assert (elem->level() == this->level());
#endif
  elem->subdomain_id() = this->subdomain_id();
  elem->processor_id() = this->processor_id();
  elem->set_id()       = this->id();
  
  // Assign the connectivity
  libmesh_assert (elem->n_nodes() == this->n_nodes());

  for (unsigned int n=0; n<elem->n_nodes(); n++)
    elem->set_node(n) = mesh.node_ptr (this->node(n));
  
  return elem;
}