dof_map.h

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

// $Id: dof_map.h 2795 2008-04-14 03:05:20Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef __dof_map_h__
#define __dof_map_h__

// C++ Includes   -----------------------------------
#include <vector>
#include <map>
#include <string>
#include <algorithm>
#include <iterator>

// Local Includes -----------------------------------
#include "libmesh_common.h"
#include "enum_order.h"
#include "reference_counted_object.h"
#include "libmesh.h" // libMesh::invalid_uint
#include "vector_value.h" // RealVectorValue
#include "threads.h"
#include "threads_allocators.h"
#include "elem_range.h"


// Forward Declarations
class DofMap;
class DofObject;
class Elem;
class MeshBase;
class Mesh;
class FEType;
class CouplingMatrix;
class System;
template <typename T> class DenseVectorBase;
template <typename T> class DenseVector;
template <typename T> class DenseMatrix;
template <typename T> class SparseMatrix;
template <typename T> class NumericVector;




// ------------------------------------------------------------
// Sparsity Pattern

/**
 * This defines the sparsity pattern, or graph, of a sparse matrix.
 * The format is quite simple -- the global indices of the nonzero entries
 * in each row are packed into a vector.  The global indices (i,j) of the 
 * nth nonzero entry of row i are given by j = sparsity_pattern[i][n];
 */
namespace SparsityPattern // use a namespace so member classes can be forward-declared.
{
  typedef std::vector<unsigned int, Threads::scalable_allocator<unsigned int> > Row;
  class Graph : public std::vector<Row> {};
  
  /**
   * Splices the two sorted ranges [begin,middle) and [middle,end)
   * into one sorted range [begin,end).  This method is much like
   * std::inplace_merge except it assumes the intersection
   * of the two sorted ranges is empty and that any element in
   * each range occurs only once in that range.  Additionally,
   * this sort occurs in-place, while std::inplace_merge may
   * use a temporary buffer.
   */ 
  template<typename BidirectionalIterator>
  static void sort_row (const BidirectionalIterator begin,
			BidirectionalIterator       middle,
			const BidirectionalIterator end);

  /**
   * This helper class can be called on multiple threads to compute 
   * the sparsity pattern (or graph) of the sparse matrix resulting
   * from the discretization.  This pattern may be used directly by
   * a particular sparse matrix format (e.g. \p LaspackMatrix)
   * or indirectly (e.g. \p PetscMatrix).  In the latter case the
   * number of nonzeros per row of the matrix is needed for efficient 
   * preallocation.  In this case it suffices to provide estimate
   * (but bounding) values, and in this case the threaded method can
   * take some short-cuts for efficiency.
   */
  class Build
  {
  private:
    const MeshBase &mesh;
    const DofMap &dof_map;
    const CouplingMatrix *dof_coupling;
    const bool implicit_neighbor_dofs;
    const bool need_full_sparsity_pattern;

  public:

    SparsityPattern::Graph sparsity_pattern;
    std::vector<unsigned int> n_nz;
    std::vector<unsigned int> n_oz;
    
    Build (const MeshBase &mesh_in,
	   const DofMap &dof_map_in,
	   const CouplingMatrix *dof_coupling_in,
	   const bool implicit_neighbor_dofs_in,
	   const bool need_full_sparsity_pattern_in) :
      mesh(mesh_in),
      dof_map(dof_map_in),
      dof_coupling(dof_coupling_in),
      implicit_neighbor_dofs(implicit_neighbor_dofs_in),
      need_full_sparsity_pattern(need_full_sparsity_pattern_in)
    {}

    Build (Build &other, Threads::split) :
      mesh(other.mesh),
      dof_map(other.dof_map),
      dof_coupling(other.dof_coupling),
      implicit_neighbor_dofs(other.implicit_neighbor_dofs),
      need_full_sparsity_pattern(other.need_full_sparsity_pattern)
    {}

    void operator()(const ConstElemRange &range);

    void join (const Build &other);
  };

#if defined(__GNUC__) && (__GNUC__ < 4) && !defined(__INTEL_COMPILER)
  /**
   * Dummy function that does nothing but can be used to prohibit
   * compiler optimization in some situations where some compilers
   * have optimization bugs.
   */
  void _dummy_function(void);
#endif
  
}




// ------------------------------------------------------------
// AMR constraint matrix types

#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)
/**
 * A row of the Dof constraint matrix.
 */
typedef std::map<unsigned int, Real, 
                 std::less<unsigned int>, 
                 Threads::scalable_allocator<std::pair<const unsigned int, Real> > > DofConstraintRow;

/** 
 * The constraint matrix storage format. 
 * We're using a class instead of a typedef to allow forward
 * declarations and future flexibility.  Don't delete this from
 * a pointer-to-std::map; the destructor isn't virtual!
 */
class DofConstraints : public std::map<unsigned int, 
                                       DofConstraintRow, 
                                       std::less<unsigned int>, 
                                       Threads::scalable_allocator<std::pair<const unsigned int, DofConstraintRow> > >
{
};
#endif // ENABLE_AMR || ENABLE_PERIODIC

  
// ------------------------------------------------------------
// Periodic boundary conditions information

#ifdef ENABLE_PERIODIC
/**
 * The definition of a periodic boundary.
 */
class PeriodicBoundary
{
public:
  // The boundary ID of this boundary and it's counterpart
  unsigned int myboundary,
	       pairedboundary;

  // One of these days we'll support rotated boundaries
  // RealTensor rotation_matrix;

  // The vector which is added to points in myboundary
  // to produce corresponding points in pairedboundary
  RealVectorValue translation_vector;
};


/** 
 * The constraint matrix storage format. 
 * We're using a class instead of a typedef to allow forward
 * declarations and future flexibility.  Is there some issue with
 * deriving from standard containers, i.e. don't do it because they
 * don't have virtual destructors?
 */
class PeriodicBoundaries : public std::map<unsigned int, PeriodicBoundary>
{
public:
  PeriodicBoundary *boundary(unsigned int id);

  PeriodicBoundaries() {}

  ~PeriodicBoundaries();

  // The periodic neighbor of \p e in direction \p side, if it
  // exists.  NULL otherwise
  const Elem *neighbor(unsigned int boundary_id, const MeshBase &mesh, const Elem *e, unsigned int side);

private:
};
#endif // ENABLE_PERIODIC


  
// ------------------------------------------------------------
// Dof Map class definition

/**
 * This class handles the numbering of degrees of freedom on a mesh.
 * For systems of equations the class supports a fixed number of variables.
 * The degrees of freedom are numbered such that sequential, contiguous blocks
 * correspond to distinct subdomains.  This is so that the resulting data
 * structures will work well with parallel linear algebra packages.
 *
 * @author Benjamin S. Kirk, 2002-2007
 */
class DofMap : public ReferenceCountedObject<DofMap>
{
public:

  /**
   * Constructor.  Requires the number of the system for which we
   * will be numbering degrees of freedom.
   */
  DofMap(const unsigned int sys_number);

  /**
   * Destructor.
   */
  ~DofMap();

  /**
   * Additional matrices may be handled with this \p DofMap.
   * They are initialized to the same sparsity structure as
   * the major matrix.
   */
  void attach_matrix (SparseMatrix<Number>& matrix);

  /**
   * Distrubute dofs on the current mesh.  Also builds the send list for
   * processor \p proc_id, which defaults to 0 for ease of use in serial
   * applications. 
   */
  void distribute_dofs (MeshBase&);

  /**
   * Computes the sparsity pattern for the matrix corresponding
   * to \p proc_id.  Produces data that can be fed to Petsc for
   * preallocation of sparse matrices.
   */
  void compute_sparsity (const MeshBase&);

  /**
   * Returns a constant reference to the \p _send_list for this processor.  The
   * \p _send_list contains the global indices of all the variables in the
   * global solution vector that influence the current processor.  This
   * information can be used for gathers at each solution step to retrieve
   * solution values needed for computation.
   */
  const std::vector<unsigned int>& get_send_list() const { return _send_list; }
  
  /**
   * Returns a constant reference to the \p _n_nz list for this processor.
   * The vector contains the bandwidth of the on-processor coupling for each
   * row of the global matrix that the current processor owns.  This
   * information can be used to preallocate space for a parallel sparse matrix.
   */
  const std::vector<unsigned int>& get_n_nz() const { return _n_nz; }
  
  /**
   * Returns a constant reference to the \p _n_oz list for this processor.
   * The vector contains the bandwidth of the off-processor coupling for each
   * row of the global matrix that the current processor owns.  This
   * information can be used to preallocate space for a parallel sparse matrix.
   */
  const std::vector<unsigned int>& get_n_oz() const { return _n_oz; }

  /**
   * Add an unknown of order \p order and finite element type
   * \p type to the system of equations.
   */
  void add_variable (const FEType& type);
  
  /**
   * @returns the approximation order for variable \p c.
   */
  Order variable_order (const unsigned int c) const;
  
  /**
   * @returns the finite element type for variable \p c.
   */
  const FEType& variable_type (const unsigned int c) const
  { return *_variable_types[c]; }
  
  /**
   * Returns the number of variables in the global solution vector. Defaults
   * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
   * Stokes (u,v,p), etc...
   */  
  unsigned int n_variables() const
  { return _variable_types.size(); }
  
  /**
   * @returns the total number of degrees of freedom in the problem.
   */
  unsigned int n_dofs() const { return _n_dfs; }

  /**
   * @returns the number of degrees of freedom on this processor.
   */
  unsigned int n_local_dofs () const
  { return this->n_dofs_on_processor (libMesh::processor_id()); }
    
  /**
   * Returns the number of degrees of freedom on subdomain \p proc.
   */
  unsigned int n_dofs_on_processor(const unsigned int proc) const
  { libmesh_assert(proc < _first_df.size()); return (_end_df[proc] - _first_df[proc]); }

  /**
   * Returns the first dof index that is local to subdomain \p proc.
   */
  unsigned int first_dof(const unsigned int proc = libMesh::processor_id()) const
  { libmesh_assert(proc < _first_df.size()); return _first_df[proc]; }
  
  /**
   * Returns the last dof index that is local to subdomain \p proc.
   * This function is now deprecated, because it returns nonsense in the rare
   * case where \p proc has no local dof indices.  Use end_dof() instead.
   */
  unsigned int last_dof(const unsigned int proc = libMesh::processor_id()) const
  { libmesh_assert(proc < _end_df.size()); return (_end_df[proc] - 1); }  

  /**
   * Returns the first dof index that is after all indices local to subdomain \p proc.
   * Analogous to the end() member function of STL containers.
   */
  unsigned int end_dof(const unsigned int proc = libMesh::processor_id()) const
  { libmesh_assert(proc < _end_df.size()); return _end_df[proc]; }  
  
  /**
   * Fills the vector \p di with the global degree of freedom indices
   * for the element. If no variable number is specified then all
   * variables are returned.
   */
  void dof_indices (const Elem* const elem,
		    std::vector<unsigned int>& di,
		    const unsigned int vn = libMesh::invalid_uint) const;

  /**
   * Builds the local element vector \p Ue from the global vector \p Ug,
   * accounting for any constrained degrees of freedom.  For an element
   * without constrained degrees of freedom this is the trivial mapping
   * \f$ Ue[i] = Ug[dof_indices[i]] \f$
   *
   * Note that the user must ensure that the element vector \p Ue is
   * properly sized when calling this method.  This is because there
   * is no \p resize() method in the \p DenseVectorBase<> class. 
   */
  void extract_local_vector (const NumericVector<Number>& Ug,
			     const std::vector<unsigned int>& dof_indices,
			     DenseVectorBase<Number>& Ue) const;


#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)

  //--------------------------------------------------------------------
  // Constraint-specific methods
  /**
   * @returns the total number of constrained degrees of freedom
   * in the problem.
   */
  unsigned int n_constrained_dofs() const { return _dof_constraints.size(); }

  /**
   * Rebuilds the raw degree of freedom constraints.
   */ 
  void create_dof_constraints (const MeshBase& mesh);

  /**
   * Gathers any relevant constraint equations from other processors
   */
  void allgather_recursive_constraints ();

  /**
   * Postprocesses any constrained degrees of freedom in elem_dofs
   * to be constrained only in terms of unconstrained dofs.
   */
  void process_recursive_constraints ();

  /**
   * Adds a copy of the user-defined row to the constraint matrix.
   * By default, produces an error if the DOF was already constrained.
   */
  void add_constraint_row (const unsigned int dof_number,
			   const DofConstraintRow& constraint_row,
			   const bool forbid_constraint_overwrite = true);
  
  /**
   * @returns true if the degree of freedom dof is constrained,
   * false otherwise.
   */
  bool is_constrained_dof (const unsigned int dof) const;
  
  /**
   * Constrains the element matrix.  This method requires the
   * element matrix to be square, in which case the elem_dofs
   * correspond to the global DOF indices of both the rows and
   * columns of the element matrix.  For this case the rows
   * and columns of the matrix necessarily correspond to variables
   * of the same approximation order.
   *
   * If \p asymmetric_constraint_rows is set to true (as it is by
   * default), constraint row equations will be reinforced in a way
   * which breaks matrix symmetry but makes inexact linear solver
   * solutions more likely to satisfy hanging node constraints.
   */
  void constrain_element_matrix (DenseMatrix<Number>& matrix,
				 std::vector<unsigned int>& elem_dofs,
				 bool asymmetric_constraint_rows = true) const;
  
  /**
   * Constrains the element matrix.  This method allows the
   * element matrix to be non-square, in which case the row_dofs
   * and col_dofs may be of different size and correspond to
   * variables approximated in different spaces.
   */
  void constrain_element_matrix (DenseMatrix<Number>& matrix,
				 std::vector<unsigned int>& row_dofs,
				 std::vector<unsigned int>& col_dofs,
				 bool asymmetric_constraint_rows = true) const;
  
  /**