coupling_matrix.h

浏览该文件的文档.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2023 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// 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 LIBMESH_COUPLING_MATRIX_H
#define LIBMESH_COUPLING_MATRIX_H

// Local Includes
#include "libmesh/libmesh_common.h"

// C++ includes
#include <algorithm>
#include <limits>
#include <utility> // std::pair
#include <vector>

namespace libMesh
{

// Forward declarations
class ConstCouplingAccessor;
class CouplingAccessor;

class ConstCouplingRow;

class ConstCouplingRowConstIterator;

class CouplingMatrix
{
public:
  explicit
  CouplingMatrix (const std::size_t n=0);

  bool operator() (const std::size_t i,
                   const std::size_t j) const;

  CouplingAccessor operator() (const std::size_t i,
                               const std::size_t j);

  std::size_t size() const;

  void resize(const std::size_t n);

  void clear();

  bool empty() const;

  CouplingMatrix & operator&= (const CouplingMatrix & other);

private:
  friend class ConstCouplingAccessor;
  friend class CouplingAccessor;
  friend class ConstCouplingRow;
  friend class ConstCouplingRowConstIterator;

  typedef std::pair<std::size_t, std::size_t> range_type;
  typedef std::vector<range_type> rc_type;
  rc_type _ranges; // 存储非零元素的范围

  std::size_t _size;
};



class ConstCouplingAccessor
{
public:
  ConstCouplingAccessor(std::size_t loc_in,
                        const CouplingMatrix & mat_in) :
    _location(loc_in), _mat(mat_in)
  {
    libmesh_assert_less(_location, _mat.size() * _mat.size());
  }

  operator bool() const
  {
    const std::size_t max_size = std::numeric_limits<std::size_t>::max();

    // 寻找可能包含i,j的范围
    CouplingMatrix::rc_type::const_iterator lb = std::upper_bound
      (_mat._ranges.begin(), _mat._ranges.end(),
       std::make_pair(_location, max_size));
    if (lb!=_mat._ranges.begin())
      --lb;
    else
      lb=_mat._ranges.end();

    // 如果没有范围可能包含i,j，则为0
    if (lb == _mat._ranges.end())
      return false;

    const std::size_t lastloc  = lb->second;

#ifdef DEBUG
    const std::size_t firstloc = lb->first;
    libmesh_assert_less_equal(firstloc, lastloc);
    libmesh_assert_less_equal(firstloc, _location);

    CouplingMatrix::rc_type::const_iterator next = lb;
    next++;
    if (next != _mat._ranges.end())
      {
        // 范围应该是排序的，不应该接触
        libmesh_assert_greater(next->first, lastloc+1);
      }
#endif

    return (lastloc >= _location);
  }

protected:
  std::size_t _location;
  const CouplingMatrix & _mat;
};


class CouplingAccessor : public ConstCouplingAccessor
{
public:
  CouplingAccessor(std::size_t loc_in,
                   CouplingMatrix & mat_in) :
    ConstCouplingAccessor(loc_in, mat_in), _my_mat(mat_in) {}

  template <typename T>
  CouplingAccessor & operator = (T new_value)
  {
    // 为了向后兼容，我们接受整数参数，但耦合矩阵条目实际上都是零或一。
    const bool as_bool = new_value;
    libmesh_assert_equal_to(new_value, as_bool);

    *this = as_bool;
    return *this;
  }

  CouplingAccessor & operator = (bool new_value)
  {
    const std::size_t max_size = std::numeric_limits<std::size_t>::max();

    // 寻找可能包含i,j的范围
    CouplingMatrix::rc_type::iterator lb =
      std::upper_bound (_my_mat._ranges.begin(), _my_mat._ranges.end(),
                        std::make_pair(_location, max_size));
    if (lb!=_my_mat._ranges.begin())
      --lb;
    else
      lb=_my_mat._ranges.end();

    // 如果没有范围可能包含i,j，则可能需要创建一个新的范围。
    if (lb == _my_mat._ranges.end())
      {
        if (new_value == true)
          _my_mat._ranges.emplace
            (_my_mat._ranges.begin(), _location, _location);
      }
    else
      {
        const std::size_t firstloc = lb->first;
        const std::size_t lastloc  = lb->second;
        libmesh_assert_less_equal(firstloc, lastloc);
        libmesh_assert_less_equal(firstloc, _location);

#ifdef DEBUG
        {
          CouplingMatrix::rc_type::const_iterator next = lb;
          next++;
          if (next != _my_mat._ranges.end())
            {
              // 范围应该是排序的，不应该接触
              libmesh_assert_greater(next->first, lastloc+1);
            }
        }
#endif

         // 如果在此范围内，可能需要缩短、移除或拆分它
    if (new_value == false)
    {
        // 如果位置在范围的起始位置，且为唯一元素，则移除整个范围
        if (_location == firstloc)
        {
            if (_location == lastloc)
            {
                _my_mat._ranges.erase(lb);
            }
            else
            {
                libmesh_assert_less(lb->first, lastloc);
                lb->first++;
            }
        }
        // 如果位置在范围的结束位置，则减小范围的结束位置
        else if (_location == lastloc)
        {
            libmesh_assert_less(firstloc, lb->second);
            lb->second--;
        }
        // 如果位置在范围内部，则调整范围的起始位置和添加新的范围
        else if (_location < lastloc)
        {
            libmesh_assert_less_equal(_location + 1, lastloc);

            lb->first = _location + 1;

            libmesh_assert_less_equal(firstloc, _location - 1);

            _my_mat._ranges.emplace(lb, firstloc, _location - 1);
        }
    }
    // 如果不在范围内，可能需要扩展范围、与相邻范围连接或添加新范围
    else // new_value == true
    {
        CouplingMatrix::rc_type::iterator next = lb;
        next++;
        const std::size_t nextloc =
            (next == _my_mat._ranges.end()) ?
            std::numeric_limits<std::size_t>::max() :
            next->first;

        // 范围应该是排序的，不应该相邻
        libmesh_assert_greater(nextloc, lastloc + 1);

        // 如果位置在范围之外，可能需要扩展范围或与下一个相邻范围连接
        if (_location > lastloc)
        {
            // 如果位置在当前范围的结束位置之后，且与下一个范围的开始位置相邻，则将其与下一个范围合并
            if (_location == lastloc + 1)
            {
                if (_location == nextloc - 1)
                {
                    next->first = firstloc;
                    _my_mat._ranges.erase(lb);
                }
                else
                    lb->second++;
            }
            // 如果位置在下一个范围的开始位置之前，则将其作为新范围插入
            else
            {
                // 如果位置刚好在下一个范围的开始位置之前，则减小下一个范围的开始位置
                if (_location == nextloc - 1)
                    next->first--;
                else
                    _my_mat._ranges.emplace(next, _location, _location);
            }
        }
    }

    return *this;
}

private:

CouplingMatrix & _my_mat;
};



class ConstCouplingRow
{
public:
    ConstCouplingRow(std::size_t row_in,
                     const CouplingMatrix & mat_in) :
        _row_i(row_in), _mat(mat_in),
        _end_location(_row_i * _mat.size() + _mat.size() - 1)
    {
        libmesh_assert_less(_row_i, _mat.size());

        // 为 i,N 查找开始范围的位置
        _begin_location = _end_location;

        const std::size_t max_size = std::numeric_limits<std::size_t>::max();

        // 查找可能包含 i,N 的范围，使用 upper_bound 不完全符合需求
        _begin_it = std::upper_bound
            (_mat._ranges.begin(), _mat._ranges.end(),
             std::make_pair(_begin_location, max_size));
        if (_begin_it != _mat._ranges.begin())
            --_begin_it;
        else
            _begin_it = _mat._ranges.end();

        // 如果范围不存在，则表示为空行
        if (_begin_it == _mat._ranges.end())
            _begin_location = max_size;
        else
        {
            const std::size_t lastloc = _begin_it->second;
            // 如果范围在 i,0 之前结束，则为空行
            std::size_t zero_location = _row_i * _mat.size();
            if (zero_location > lastloc)
            {
                _begin_location = max_size;
                _begin_it = _mat._ranges.end();
            }
            else
            {
                while (_begin_it != _mat._ranges.begin())
                {
                    CouplingMatrix::rc_type::const_iterator prev =
                        _begin_it;
                    --prev;

                    if (prev->second < zero_location)
                        break;

                    _begin_it = prev;
                }
                if (_begin_it->first < zero_location)
                    _begin_location = zero_location;
                else
                    _begin_location = _begin_it->first;
            }
        }
    }

    /*
    * \brief 用于在此行中循环遍历索引的前向迭代器类型。
    */
    typedef ConstCouplingRowConstIterator const_iterator;

    /*
     * \brief 返回此行中第一个索引的迭代器，对于空行则返回 end()。
     */
    const_iterator begin() const;

    /*
     * \brief 表示行末尾的迭代器。
     */
    const_iterator end() const;

    bool operator== (const ConstCouplingRow & other) const
    {
        // 考虑到来自不同矩阵对象的行是相等的是不正确的
        libmesh_assert(&_mat == &other._mat);

        return ((_begin_location == other._begin_location) &&
            (_begin_it == other._begin_it));
    }

    bool operator!= (const ConstCouplingRow & other) const
    {
        return !(*this == other);
    }
protected:

    friend class ConstCouplingRowConstIterator;

    std::size_t _row_i;
    const CouplingMatrix & _mat;

    // 对应于此行第一个条目的位置 (i*size+j)，对于空行则为 std::numeric_limits<size_t>::max()
    std::size_t _begin_location;

    // 对应于此行结束的位置 (i*size+j)（不考虑该结束位置是否在范围内）。缓存位置，因为每次 ++iterator 调用都需要使用它
    const std::size_t _end_location;

    // 包含第一个行元素的范围的迭代器，如果对于此行不存在 CouplingMatrix 值，则为 _mat._ranges.end()
    CouplingMatrix::rc_type::const_iterator _begin_it;
};



class ConstCouplingRowConstIterator
{
public:
    ConstCouplingRowConstIterator(const ConstCouplingRow & row_in,
        std::size_t loc_in,
        CouplingMatrix::rc_type::const_iterator it_in) :
        _location(loc_in),
        _row(row_in),
        _it(it_in)
    {
#ifndef NDEBUG
        if (_it != _row._mat._ranges.end())
        {
            libmesh_assert_less_equal(_it->first, _location);
            libmesh_assert_less_equal(_location, _it->second);
        }
        else
        {
            libmesh_assert_equal_to
                (_location, std::numeric_limits<size_t>::max());
        }
#endif
    }

    std::size_t operator* ()
    {
        libmesh_assert_not_equal_to
            (_location, std::numeric_limits<std::size_t>::max());
        return _location % _row._mat.size();
    }

    ConstCouplingRowConstIterator & operator++ ()
    {
        libmesh_assert_not_equal_to
            (_location, std::numeric_limits<std::size_t>::max());
        libmesh_assert
            (_it != _row._mat._ranges.end());

        if (_location >= _it->second)
        {
            ++_it;

            // 是否超出矩阵的末尾？
            if (_it == _row._mat._ranges.end())
                _location = std::numeric_limits<std::size_t>::max();
            else
            {
                _location = _it->first;
                // 新范围是否超出行的末尾？
                if (_location > _row._end_location)
                    _location = std::numeric_limits<std::size_t>::max();
            }
        }
        // 递增后是否超出行的末尾？
        else if (_location >= _row._end_location)
            _location = std::numeric_limits<std::size_t>::max();
        else
            ++_location;

        return *this;
    }

    bool operator== (const ConstCouplingRowConstIterator & other) const
    {
        // 考虑到来自不同行对象的迭代器是相等的是不正确的
        libmesh_assert(_row == other._row);

        return (_location == other._location);
    }

    bool operator!= (const ConstCouplingRowConstIterator & other) const
    {
        return !(*this == other);
    }

private:
    // 对应于此迭代器位置的位置 (i*size+j)，或 std::numeric_limits<size_t>::max() 表示 end()
    std::size_t _location;
    const ConstCouplingRow & _row;
    // 包含此迭代器位置的范围，或 _row._mat._ranges.end() 表示此位置不在任何范围内
    CouplingMatrix::rc_type::const_iterator _it;
};



//--------------------------------------------------
// ConstCouplingRow inline methods
inline
ConstCouplingRow::const_iterator ConstCouplingRow::begin() const
{
  return const_iterator (*this, _begin_location, _begin_it);
}

inline
ConstCouplingRow::const_iterator ConstCouplingRow::end() const
{
  return const_iterator
    (*this, std::numeric_limits<std::size_t>::max(),
     _mat._ranges.end());
}



//--------------------------------------------------
// CouplingMatrix inline methods
inline
CouplingMatrix::CouplingMatrix (const std::size_t n) :
  _ranges(), _size(n)
{
  this->resize(n);
}



inline
bool CouplingMatrix::operator() (const std::size_t i,
                                 const std::size_t j) const
{
  libmesh_assert_less (i, _size);
  libmesh_assert_less (j, _size);

  const std::size_t location = std::size_t(i)*_size + j;

  return bool(ConstCouplingAccessor(location, *this));
}




inline
CouplingAccessor CouplingMatrix::operator() (const std::size_t i,
                                             const std::size_t j)
{
  const std::size_t location = std::size_t(i)*_size + j;

  return CouplingAccessor(location, *this);
}



inline
std::size_t CouplingMatrix::size() const
{
  return _size;
}



inline
void CouplingMatrix::resize(const std::size_t n)
{
  _size = n;

  _ranges.clear();
}



inline
void CouplingMatrix::clear()
{
  this->resize(0);
}



inline
bool CouplingMatrix::empty() const
{
  return (_size == 0);
}


} // namespace libMesh


#endif // LIBMESH_COUPLING_MATRIX_H