Loading...
Searching...
No Matches
form_vf_assembly.cc
Go to the documentation of this file.
//
#include "rheolef/form_vf_assembly.h"
#include "rheolef/eigen_util.h"
namespace rheolef { namespace details {
using namespace std;
// ------------------------------------------------------------------------------
// norm_max
// ------------------------------------------------------------------------------
template <class T>
T
norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
T m_max = 0;
for (size_t i = 0; i < size_t(m.rows()); i++) {
for (size_t j = 0; j < size_t(m.cols()); j++) {
m_max = std::max (m_max, m(i,j));
}}
return m_max;
}
// ------------------------------------------------------------------------------
// is_symmetric
// ------------------------------------------------------------------------------
template <class T>
bool
check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max) {
if (m.rows() != m.cols()) return false; // non-square matrix
const T eps = std::numeric_limits<T>::epsilon();
if (tol_m_max < sqr(eps)) return true; // zero matrix
for (size_t i = 0; i < size_t(m.rows()); i++) {
for (size_t j = i+1; j < size_t(m.cols()); j++) {
if (abs(m(i,j) - m(j,i)) > tol_m_max) return false;
}}
return true;
}
// ------------------------------------------------------------------------------
// local lump
// ------------------------------------------------------------------------------
template <class T>
void
local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
check_macro (m.rows() == m.cols(), "unexpected rectangular matrix for lumped mass");
for (size_t i = 0; i < size_t(m.rows()); i++) {
T sum = 0;
for (size_t j = 0; j < size_t(m.cols()); j++) {
sum += m(i,j);
m(i,j) = T(0);
}
m(i,i) = sum;
}
}
// ------------------------------------------------------------------------------
// local inversion
// ------------------------------------------------------------------------------
template <class T>
void
local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, bool is_diag) {
check_macro (m.rows() == m.cols(), "unexpected rectangular matrix for local invert");
if (is_diag) {
// matrix has been lumped: easy diagonal inversion
for (size_t i = 0; i < size_t(m.rows()); i++) {
m(i,i) = 1./m(i,i);
}
return;
}
// general inversion
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> inv_m (m.rows(),m.cols());
bool status = invert(m,inv_m);
if (!status) {
#ifndef TO_CLEAN
std::ofstream out ("inv_failed.mtx");
put_matrix_market (out, m);
out.close();
#endif // TO_CLEAN
error_macro("singular element matrix founded");
}
m = inv_m;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciate(T) \
template T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m); \
template bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, const T&); \
template void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&); \
template void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, bool); \
_RHEOLEF_instanciate(Float)
}} // namespace rheolef::details
Float
see the Float page for the full documentation
_RHEOLEF_instanciate
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::details
Definition compiler_eigen.h:39
rheolef::details::local_lump
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition form_vf_assembly.cc:61
rheolef::details::check_is_symmetric
bool check_is_symmetric(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, const T &tol_m_max)
Definition form_vf_assembly.cc:46
rheolef::details::norm_max
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition form_vf_assembly.cc:33
rheolef::details::local_invert
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
Definition form_vf_assembly.cc:77
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
rheolef::put_matrix_market
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition eigen_util.h:78