dof_map_constraints.c

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

// $Id: dof_map_constraints.C 2789 2008-04-13 02:24:40Z 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



// C++ Includes -------------------------------------
#include <set>
#include <algorithm> // for std::count, std::fill

// Local Includes -----------------------------------
#include "dof_map.h"
#include "elem.h"
#include "mesh_base.h"
#include "fe_interface.h"
#include "fe_type.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "libmesh_logging.h"
#include "system.h" // needed by enforce_constraints_exactly()
#include "mesh.h"   // as is this
#include "numeric_vector.h" // likewise
#include "parallel.h"
#include "point_locator_base.h"
#include "elem_range.h"
#include "threads_allocators.h"



// Anonymous namespace to hold helper classes
namespace {

  class ComputeConstraints
  {
  public:
    ComputeConstraints (DofConstraints &constraints,
			DofMap &dof_map,
#ifdef ENABLE_PERIODIC
			PeriodicBoundaries &periodic_boundaries,
#endif
			const MeshBase &mesh,
			const unsigned int variable_number) :
      _constraints(constraints),
      _dof_map(dof_map),
#ifdef ENABLE_PERIODIC
      _periodic_boundaries(periodic_boundaries),
#endif
      _mesh(mesh),
      _variable_number(variable_number)
    {}

    void operator()(const ConstElemRange &range) const
    {
      for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
        {
#ifdef ENABLE_AMR
	  FEInterface::compute_constraints (_constraints,
					    _dof_map,
					    _variable_number,
					    *it);
#endif
#ifdef ENABLE_PERIODIC
          // FIXME: periodic constraints won't work on a non-serial
          // mesh unless it's kept ghost elements from opposing
          // boundaries!
	  FEInterface::compute_periodic_constraints (_constraints,
					             _dof_map,
                                                     _periodic_boundaries,
						     _mesh,
					             _variable_number,
					             *it);
#endif
        }
    }

  private:
    DofConstraints &_constraints;
    DofMap &_dof_map;
#ifdef ENABLE_PERIODIC
    PeriodicBoundaries &_periodic_boundaries;
#endif
    const MeshBase &_mesh;
    const unsigned int _variable_number;    
  };
}



// ------------------------------------------------------------
// DofMap member functions

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

void DofMap::create_dof_constraints(const MeshBase& mesh)
{
  START_LOG("create_dof_constraints()", "DofMap");

  libmesh_assert (mesh.is_prepared());
  
  const unsigned int dim = mesh.mesh_dimension();

  // Constraints are not necessary in 1D
  if (dim == 1)
  {
    // make sure we stop logging though
    STOP_LOG("create_dof_constraints()", "DofMap");
    return;
  }
  
  // Here we build the hanging node constraints.  This is done
  // by enforcing the condition u_a = u_b along hanging sides.
  // u_a = u_b is collocated at the nodes of side a, which gives
  // one row of the constraint matrix.
  
  // clear any existing constraints.
  _dof_constraints.clear();

  // define the range of elements of interest
  ConstElemRange range;
  {
    // With SerialMesh or a serial ParallelMesh, every processor
    // computes every constraint
    MeshBase::const_element_iterator
      elem_begin = mesh.elements_begin(),
      elem_end   = mesh.elements_end();
    
    // With a parallel ParallelMesh, processors compute only
    // their local constraints
    if (!mesh.is_serial())
      {
	elem_begin = mesh.local_elements_begin();
	elem_end   = mesh.local_elements_end(); 
      }

    // set the range to contain the specified elements
    range.reset (elem_begin, elem_end);
  }
  
  // Look at all the variables in the system.  Reset the element
  // range at each iteration -- there is no need to reconstruct it.
  for (unsigned int variable_number=0; variable_number<this->n_variables();
       ++variable_number, range.reset())
    Threads::parallel_for (range,
			   ComputeConstraints (_dof_constraints,
					       *this,
#ifdef ENABLE_PERIODIC
					       _periodic_boundaries,
#endif
					       mesh,
					       variable_number));

  // With a parallelized Mesh, we've computed our local constraints,
  // but they may depend on non-local constraints that we'll need to
  // take into account.
  if (!mesh.is_serial())
    this->allgather_recursive_constraints();
  
  STOP_LOG("create_dof_constraints()", "DofMap");
}



void DofMap::add_constraint_row (const unsigned int dof_number,
				 const DofConstraintRow& constraint_row,
				 const bool forbid_constraint_overwrite)
{
  // Optionally allow the user to overwrite constraints.  Defaults to false.
  if (forbid_constraint_overwrite)
    if (this->is_constrained_dof(dof_number))
      {
	std::cerr << "ERROR: DOF " << dof_number << " was already constrained!"
		  << std::endl;
	libmesh_error();
      }
  
  std::pair<unsigned int, DofConstraintRow> kv(dof_number, constraint_row);

  _dof_constraints.insert(kv);
}



void DofMap::print_dof_constraints(std::ostream& os) const
{
  os << "DOF CONSTRAINTS OUTPUT:"
     << std::endl;
  
  for (DofConstraints::const_iterator it=_dof_constraints.begin();
       it != _dof_constraints.end(); ++it)
    {
      const unsigned int i = it->first;
      const DofConstraintRow& row = it->second;

      os << "Constraints for DOF " << i
	 << ": \t";

      for (DofConstraintRow::const_iterator pos=row.begin();
	   pos != row.end(); ++pos)
	os << " (" << pos->first << ","
	   << pos->second << ")\t";

      os << std::endl;
    }
}



void DofMap::constrain_element_matrix (DenseMatrix<Number>& matrix,
				       std::vector<unsigned int>& elem_dofs,
				       bool asymmetric_constraint_rows) const
{
  libmesh_assert (elem_dofs.size() == matrix.m());
  libmesh_assert (elem_dofs.size() == matrix.n());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
    
  // The constrained matrix is built up as C^T K C.    
  DenseMatrix<Number> C;

  
  this->build_constraint_matrix (C, elem_dofs);

  START_LOG("constrain_elem_matrix()", "DofMap");
  
  // It is possible that the matrix is not constrained at all.
  if ((C.m() == matrix.m()) &&
      (C.n() == elem_dofs.size())) // It the matrix is constrained
    {
      // Compute the matrix-matrix-matrix product C^T K C
      matrix.left_multiply_transpose  (C);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == matrix.n());
      libmesh_assert (matrix.m() == elem_dofs.size());
      libmesh_assert (matrix.n() == elem_dofs.size());
      
      
      for (unsigned int i=0; i<elem_dofs.size(); i++)
	// If the DOF is constrained
	if (this->is_constrained_dof(elem_dofs[i]))
	  {
	    for (unsigned int j=0; j<matrix.n(); j++)
	      matrix(i,j) = 0.;
	    
	    matrix(i,i) = 1.;
	    
            if (asymmetric_constraint_rows)
              {
	        DofConstraints::const_iterator
	          pos = _dof_constraints.find(elem_dofs[i]);
	    
	        libmesh_assert (pos != _dof_constraints.end());
	    
	        const DofConstraintRow& constraint_row = pos->second;
	    
	        libmesh_assert (!constraint_row.empty());
	    
	        for (DofConstraintRow::const_iterator
		       it=constraint_row.begin(); it != constraint_row.end();
		     ++it)
	          for (unsigned int j=0; j<elem_dofs.size(); j++)
		    if (elem_dofs[j] == it->first)
		      matrix(i,j) = -it->second;	
	      }
	  }
    } // end if is constrained...
  
  STOP_LOG("constrain_elem_matrix()", "DofMap");  
}



void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number>& matrix,
						  DenseVector<Number>& rhs,
						  std::vector<unsigned int>& elem_dofs,
						  bool asymmetric_constraint_rows) const
{
  libmesh_assert (elem_dofs.size() == matrix.m());
  libmesh_assert (elem_dofs.size() == matrix.n());
  libmesh_assert (elem_dofs.size() == rhs.size());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained matrix is built up as C^T K C.
  // The constrained RHS is built up as C^T F
  DenseMatrix<Number> C;
  
  this->build_constraint_matrix (C, elem_dofs);
  
  START_LOG("cnstrn_elem_mat_vec()", "DofMap");
  
  // It is possible that the matrix is not constrained at all.
  if ((C.m() == matrix.m()) &&
      (C.n() == elem_dofs.size())) // It the matrix is constrained
    {
      // Compute the matrix-matrix-matrix product C^T K C
      matrix.left_multiply_transpose  (C);
      matrix.right_multiply (C);
      
      
      libmesh_assert (matrix.m() == matrix.n());
      libmesh_assert (matrix.m() == elem_dofs.size());
      libmesh_assert (matrix.n() == elem_dofs.size());
      

      for (unsigned int i=0; i<elem_dofs.size(); i++)
	if (this->is_constrained_dof(elem_dofs[i]))
	  {
	    for (unsigned int j=0; j<matrix.n(); j++)
	      matrix(i,j) = 0.;
	    
	    // If the DOF is constrained
	    matrix(i,i) = 1.;
	    
            // This will put a nonsymmetric entry in the constraint
            // row to ensure that the linear system produces the
            // correct value for the constrained DOF.
            if (asymmetric_constraint_rows)
              {
	        DofConstraints::const_iterator
	          pos = _dof_constraints.find(elem_dofs[i]);
	    
	        libmesh_assert (pos != _dof_constraints.end());
	    
	        const DofConstraintRow& constraint_row = pos->second;
	    
// p refinement creates empty constraint rows
//	    libmesh_assert (!constraint_row.empty());
	    
	        for (DofConstraintRow::const_iterator
		       it=constraint_row.begin(); it != constraint_row.end();
		     ++it)
	          for (unsigned int j=0; j<elem_dofs.size(); j++)
		    if (elem_dofs[j] == it->first)
		      matrix(i,j) = -it->second;	
              }
	  }

      
      // Compute the matrix-vector product C^T F
      DenseVector<Number> old_rhs(rhs);

      // resize the RHS vector & 0 before summation
      rhs.resize(elem_dofs.size());
      rhs.zero();
      
      // compute matrix/vector product
      for (unsigned int i=0; i<elem_dofs.size(); i++)
	for (unsigned int j=0; j<old_rhs.size(); j++)
	  rhs(i) += C.transpose(i,j)*old_rhs(j);
	

      libmesh_assert (elem_dofs.size() == rhs.size());

      for (unsigned int i=0; i<elem_dofs.size(); i++)
	if (this->is_constrained_dof(elem_dofs[i]))
	  {	
	    // If the DOF is constrained
	    rhs(i) = 0.;
	  }
    } // end if is constrained...
  
  STOP_LOG("cnstrn_elem_mat_vec()", "DofMap");  
}



void DofMap::constrain_element_matrix (DenseMatrix<Number>& matrix,
				       std::vector<unsigned int>& row_dofs,
				       std::vector<unsigned int>& col_dofs,
				       bool asymmetric_constraint_rows) const
{
  libmesh_assert (row_dofs.size() == matrix.m());
  libmesh_assert (col_dofs.size() == matrix.n());

  // check for easy return
  if (this->n_constrained_dofs() == 0)
    return;
  
  // The constrained matrix is built up as R^T K C.
  DenseMatrix<Number> R;
  DenseMatrix<Number> C;

  // Safeguard against the user passing us the same
  // object for row_dofs and col_dofs.  If that is done
  // the calls to build_matrix would fail
  std::vector<unsigned int> orig_row_dofs(row_dofs);
  std::vector<unsigned int> orig_col_dofs(col_dofs);
  
  this->build_constraint_matrix (R, orig_row_dofs);