dense_matrix_base_impl.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_DENSE_MATRIX_BASE_IMPL_H
#define LIBMESH_DENSE_MATRIX_BASE_IMPL_H

// Local Includes
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector_base.h"
#include "libmesh/dense_vector.h"
#include "libmesh/int_range.h"

// C++ includes
#include <iomanip> // for std::setw()

namespace libMesh
{

template <typename T>
DenseVector<T> DenseMatrixBase<T>::diagonal() const
{
    DenseVector<T> ret(_m);
    for (decltype(_m) i = 0; i < _m; ++i)
        ret(i) = el(i, i);
    return ret;
}

template<typename T>
void DenseMatrixBase<T>::multiply(DenseMatrixBase<T> & M1,
                                   const DenseMatrixBase<T> & M2,
                                   const DenseMatrixBase<T> & M3)
{
    // 断言确保我们传递了正确维度的矩阵
    libmesh_assert_equal_to(M1.m(), M2.m());
    libmesh_assert_equal_to(M1.n(), M3.n());
    libmesh_assert_equal_to(M2.n(), M3.m());

    const unsigned int m_s = M2.m();
    const unsigned int p_s = M2.n();
    const unsigned int n_s = M1.n();

    // 通过这种方式执行是因为有很大的机会（至少对于约束矩阵来说）
    // 在右乘时 M3(k,j) = 0.
    for (unsigned int k = 0; k < p_s; k++)
        for (unsigned int j = 0; j < n_s; j++)
            if (M3.el(k, j) != 0.)
                for (unsigned int i = 0; i < m_s; i++)
                    M1.el(i, j) += M2.el(i, k) * M3.el(k, j);
}

template<typename T>
void DenseMatrixBase<T>::condense(const unsigned int iv,
                                  const unsigned int jv,
                                  const T val,
                                  DenseVectorBase<T> & rhs)
{
    libmesh_assert_equal_to(this->_m, rhs.size());
    libmesh_assert_equal_to(iv, jv);

    // 将已知值移入 RHS，并将列置零
    for (auto i : make_range(this->m()))
    {
        rhs.el(i) -= this->el(i, jv) * val;
        this->el(i, jv) = 0.;
    }

    // 将行置零
    for (auto j : make_range(this->n()))
        this->el(iv, j) = 0.;

    this->el(iv, jv) = 1.;
    rhs.el(iv) = val;
}

template<typename T>
void DenseMatrixBase<T>::print_scientific(std::ostream & os, unsigned precision) const
{
    // 保存初始格式标志
    std::ios_base::fmtflags os_flags = os.flags();

    // 打印矩阵元素
    for (auto i : make_range(this->m()))
    {
        for (auto j : make_range(this->n()))
            os << std::setw(15)
               << std::scientific
               << std::setprecision(precision)
               << this->el(i, j) << " ";

        os << std::endl;
    }

    // 恢复原始格式标志
    os.flags(os_flags);
}

template<typename T>
void DenseMatrixBase<T>::print(std::ostream & os) const
{
    for (auto i : make_range(this->m()))
    {
        for (auto j : make_range(this->n()))
            os << std::setw(8)
               << this->el(i, j) << " ";

        os << std::endl;
    }

    return;
}



}

#endif // LIBMESH_DENSE_MATRIX_BASE_IMPL_H