elem.c

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

                        if (this->contains_vertex_of(*child_it) ||
                            (*child_it)->contains_vertex_of(this))
                          neighbor_set.insert (*child_it);
                    }
#endif // #ifdef ENABLE_AMR
                }
            }
        }
    }
  while (old_size != neighbor_set.size());
}



#ifdef DEBUG

void Elem::libmesh_assert_valid_node_pointers() const
{
  libmesh_assert(this->valid_id());
  for (unsigned int n=0; n != this->n_nodes(); ++n)
    {
      libmesh_assert(this->get_node(n));
      libmesh_assert(this->get_node(n)->valid_id());
    }
}



void Elem::libmesh_assert_valid_neighbors() const
{
  for (unsigned int s=0; s<this->n_neighbors(); s++)
    {
      const Elem *neigh = this->neighbor(s);

      // Any element might have a remote neighbor; checking
      // to make sure that's not inaccurate is tough.
      if (neigh == remote_elem)
        continue;

      if (neigh)
        {
          // Only subactive elements have subactive neighbors
          libmesh_assert (this->subactive() || !neigh->subactive());

          const Elem *elem = this;

          // If we're subactive but our neighbor isn't, its
          // return neighbor link will be to our first active
          // ancestor
          if (this->subactive() && !neigh->subactive())
            {
              for (elem = this; !elem->active();
                   elem = elem->parent())
                libmesh_assert(elem);
            }

          unsigned int rev = neigh->which_neighbor_am_i(elem);
          libmesh_assert (rev < neigh->n_neighbors());

          if (this->subactive() && !neigh->subactive())
            {
              libmesh_assert (neigh->neighbor(rev) == elem);
            }
          else
            {
              Elem *nn = neigh->neighbor(rev);
              libmesh_assert(nn);

              for (; elem != nn; elem = elem->parent())
                libmesh_assert(elem);
            }
        }
      else
        {
          const Elem *parent = this->parent();
          if (parent)
            libmesh_assert (!parent->neighbor(s));
        }
    }
}

#endif // DEBUG



void Elem::make_links_to_me_remote()
{
  libmesh_assert (this != remote_elem);

  // We need to have handled any children first
#if defined(ENABLE_AMR) && defined(DEBUG)
  if (this->has_children())
    for (unsigned int c = 0; c != this->n_children(); ++c)
      {
        Elem *child = this->child(c);
        libmesh_assert (child == remote_elem);
      }
#endif

  // Remotify any neighbor links
  for (unsigned int s = 0; s != this->n_sides(); ++s)
    {
      Elem *neigh = this->neighbor(s);
      if (neigh && neigh != remote_elem)
        {
	  // My neighbor should never be more refined than me; my real
	  // neighbor would have been its parent in that case.
	  libmesh_assert(this->level() >= neigh->level());
	  
          if (this->level() == neigh->level() &&
              neigh->has_neighbor(this))
            {
#ifdef ENABLE_AMR
	      // My neighbor may have descendants which also consider me a
	      // neighbor
              std::vector<const Elem*> family;
              neigh->family_tree_by_neighbor (family, this);

              // FIXME - There's a lot of ugly const_casts here; we
              // may want to make remote_elem non-const and create
              // non-const versions of the family_tree methods
              for (unsigned int i=0; i != family.size(); ++i)
                {
                  Elem *n = const_cast<Elem*>(family[i]);
                  libmesh_assert (n);
                  if (n == remote_elem)
                    continue;
                  unsigned int my_s = n->which_neighbor_am_i(this);
                  libmesh_assert (my_s < n->n_neighbors());
                  libmesh_assert (n->neighbor(my_s) == this);
                  n->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
                }
#else
              unsigned int my_s = neigh->which_neighbor_am_i(this);
              libmesh_assert (my_s < neigh->n_neighbors());
              libmesh_assert (neigh->neighbor(my_s) == this);
              neigh->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
#endif
            }
#ifdef ENABLE_AMR
          // Even if my neighbor doesn't link back to me, it might
	  // have subactive descendants which do
	  else if (neigh->has_children())
            {
              // If my neighbor at the same level doesn't have me as a
	      // neighbor, I must be subactive
	      libmesh_assert(this->level() > neigh->level() ||
			     this->subactive());

              // My neighbor must have some ancestor of mine as a
	      // neighbor
	      Elem *ancestor = this->parent();
	      libmesh_assert(ancestor);
              while (!neigh->has_neighbor(ancestor))
                {
                  ancestor = ancestor->parent();
	          libmesh_assert(ancestor);
                }

	      // My neighbor may have descendants which consider me a
	      // neighbor
              std::vector<const Elem*> family;
              neigh->family_tree_by_subneighbor (family, ancestor, this);

              // FIXME - There's a lot of ugly const_casts here; we
              // may want to make remote_elem non-const and create
              // non-const versions of the family_tree methods
              for (unsigned int i=0; i != family.size(); ++i)
                {
                  Elem *n = const_cast<Elem*>(family[i]);
                  libmesh_assert (n);
                  if (n == remote_elem)
                    continue;
                  unsigned int my_s = n->which_neighbor_am_i(this);
                  libmesh_assert (my_s < n->n_neighbors());
                  libmesh_assert (n->neighbor(my_s) == this);
                  n->set_neighbor(my_s, const_cast<RemoteElem*>(remote_elem));
                }
            }
#endif
        }
    }

#ifdef ENABLE_AMR
  // Remotify parent's child link
  Elem *parent = this->parent();
  if (parent && parent != remote_elem)
    {
      unsigned int me = parent->which_child_am_i(this);
      libmesh_assert (parent->_children[me] == this);
      parent->_children[me] = const_cast<RemoteElem*>(remote_elem);
    }
#endif
}



void Elem::write_connectivity (std::ostream& out,
			       const IOPackage iop) const
{
  libmesh_assert (out.good());
  libmesh_assert (_nodes != NULL);
  libmesh_assert (iop != INVALID_IO_PACKAGE);

  switch (iop)
    {
    case TECPLOT:
      {
	// This connectivity vector will be used repeatedly instead
	// of being reconstructed inside the loop.
	std::vector<unsigned int> conn;
	for (unsigned int sc=0; sc <this->n_sub_elem(); sc++)
	  {
	    this->connectivity(sc, TECPLOT, conn);
	    
	    std::copy(conn.begin(),
		      conn.end(),
		      std::ostream_iterator<unsigned int>(out, " "));
	    
	    out << '\n';
	  }
	return;
      }

    case UCD:
      {
	for (unsigned int i=0; i<this->n_nodes(); i++)
	  out << this->node(i)+1 << "\t";
	
	out << '\n';
	return;
      }

    default:
      libmesh_error();
    }

  libmesh_error();
}


// void Elem::write_tecplot_connectivity(std::ostream& out) const
// {
//   libmesh_assert (!out.bad());
//   libmesh_assert (_nodes != NULL);

//   // This connectivity vector will be used repeatedly instead
//   // of being reconstructed inside the loop.
//   std::vector<unsigned int> conn;
//   for (unsigned int sc=0; sc <this->n_sub_elem(); sc++)
//     {
//       this->connectivity(sc, TECPLOT, conn);
      
//       std::copy(conn.begin(),
//  		conn.end(),
//  		std::ostream_iterator<unsigned int>(out, " "));
      
//       out << std::endl;
//     }
// }



// void Elem::write_ucd_connectivity(std::ostream &out) const
// {
//   libmesh_assert (out);
//   libmesh_assert (_nodes != NULL);

//   for (unsigned int i=0; i<this->n_nodes(); i++)
//     out << this->node(i)+1 << "\t";

//   out << std::endl;
// }



Real Elem::quality (const ElemQuality q) const
{
  switch (q)
    {    
      /**
       * I don't know what to do for this metric. 
       */
    default:
      {
	here();

	std::cerr << "ERROR:  unknown quality metric: "
		  << q 
		  << std::endl
		  << "Cowardly returning 1."
		  << std::endl;

	return 1.;
      }
    }

    
    // Will never get here...
    libmesh_error();
    return 0.;
}



bool Elem::ancestor() const
{
#ifdef ENABLE_AMR

  if (this->active())
    return false;

if (!this->has_children())
    return false;
  if (this->child(0)->active())
    return true;

  return this->child(0)->ancestor();
#else
  return false;
#endif
}



#ifdef ENABLE_AMR

void Elem::add_child (Elem* elem)
{
  if(_children == NULL)
    {
      _children = new Elem*[this->n_children()];
      
      for (unsigned int c=0; c<this->n_children(); c++)
	_children[c] = NULL;
    }
  
  for (unsigned int c=0; c<this->n_children(); c++)
    {
      if(_children[c] == NULL || _children[c] == remote_elem)
	{
	  libmesh_assert (this == elem->parent());
	  _children[c] = elem;
	  return;
	}
    }

  std::cerr << "Error: Tried to add a child to an element with full children array"
            << std::endl;
  libmesh_error();
}



void Elem::add_child (Elem* elem, unsigned int c)
{
  if(_children == NULL)
    {
      _children = new Elem*[this->n_children()];
      
      for (unsigned int i=0; i<this->n_children(); i++)
	_children[i] = NULL;
    }
  
  libmesh_assert (_children[c] == NULL || _children[c] == remote_elem);
  libmesh_assert (this == elem->parent());

  _children[c] = elem;
}



bool Elem::is_child_on_edge(const unsigned int c,
                            const unsigned int e) const
{
  libmesh_assert (c < this->n_children());
  libmesh_assert (e < this->n_edges());

  AutoPtr<Elem> my_edge = this->build_edge(e);
  AutoPtr<Elem> child_edge = this->build_edge(e);

  // We're assuming that an overlapping child edge has the same
  // number and orientation as its parent
  return (child_edge->node(0) == my_edge->node(0) ||
      child_edge->node(1) == my_edge->node(1));
}


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

  // Add this element to the family tree.
  family.push_back(this);

  // Recurse into the elements children, if it has them.
  // Do not clear the vector any more.
  if (!this->active())
    for (unsigned int c=0; c<this->n_children(); c++)
      if (!this->child(c)->is_remote())
	this->child(c)->family_tree (family, false);
}



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

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

  // Otherwise recurse into the element's children.
  // Do not clear the vector any more.
  else 
    for (unsigned int c=0; c<this->n_children(); c++)
      if (!this->child(c)->is_remote())
	this->child(c)->active_family_tree (active_family, false);
}



void Elem::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
  libmesh_assert (this->has_neighbor(neighbor));

  // Add this element to the family tree.
  family.push_back(this);

  // Recurse into the elements children, if it has them.
  // Do not clear the vector any more.
  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->family_tree_by_neighbor (family, neighbor, false);
      }
}



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

  // To simplifly this function we need an existing neighbor
  libmesh_assert (neighbor);
  libmesh_assert (neighbor != remote_elem);
  libmesh_assert (this->has_neighbor(neighbor));

  // This only makes sense if subneighbor descends from neighbor
  libmesh_assert (subneighbor);
  libmesh_assert (subneighbor != remote_elem);
  libmesh_assert (neighbor->is_ancestor_of(subneighbor));

  // Add this element to the family tree if applicable.
  if (neighbor == subneighbor)
    family.push_back(this);

  // Recurse into the elements children, if it has them.
  // Do not clear the vector any more.
  if (this->has_children())
    for (unsigned int c=0; c != this->n_children(); ++c)
      {
        Elem *child = this->child(c);
        if (child != remote_elem)