Rheolef  7.2
an efficient C++ finite element environment
form_vf_assembly.cc
Go to the documentation of this file.
1 //
22 #include "rheolef/form_vf_assembly.h"
23 #include "rheolef/eigen_util.h"
24 
25 namespace rheolef { namespace details {
26 using namespace std;
27 
28 // ------------------------------------------------------------------------------
29 // norm_max
30 // ------------------------------------------------------------------------------
31 template <class T>
32 T
33 norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
34  T m_max = 0;
35  for (size_t i = 0; i < size_t(m.rows()); i++) {
36  for (size_t j = 0; j < size_t(m.cols()); j++) {
37  m_max = std::max (m_max, m(i,j));
38  }}
39  return m_max;
40 }
41 // ------------------------------------------------------------------------------
42 // is_symmetric
43 // ------------------------------------------------------------------------------
44 template <class T>
45 bool
46 check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, const T& tol_m_max) {
47  if (m.rows() != m.cols()) return false; // non-square matrix
48  const T eps = std::numeric_limits<T>::epsilon();
49  if (tol_m_max < sqr(eps)) return true; // zero matrix
50  for (size_t i = 0; i < size_t(m.rows()); i++) {
51  for (size_t j = i+1; j < size_t(m.cols()); j++) {
52  if (abs(m(i,j) - m(j,i)) > tol_m_max) return false;
53  }}
54  return true;
55 }
56 // ------------------------------------------------------------------------------
57 // local lump
58 // ------------------------------------------------------------------------------
59 template <class T>
60 void
61 local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
62  check_macro (m.rows() == m.cols(), "unexpected rectangular matrix for lumped mass");
63  for (size_t i = 0; i < size_t(m.rows()); i++) {
64  T sum = 0;
65  for (size_t j = 0; j < size_t(m.cols()); j++) {
66  sum += m(i,j);
67  m(i,j) = T(0);
68  }
69  m(i,i) = sum;
70  }
71 }
72 // ------------------------------------------------------------------------------
73 // local inversion
74 // ------------------------------------------------------------------------------
75 template <class T>
76 void
77 local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m, bool is_diag) {
78  check_macro (m.rows() == m.cols(), "unexpected rectangular matrix for local invert");
79  if (is_diag) {
80  // matrix has been lumped: easy diagonal inversion
81  for (size_t i = 0; i < size_t(m.rows()); i++) {
82  m(i,i) = 1./m(i,i);
83  }
84  return;
85  }
86  // general inversion
87  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> inv_m (m.rows(),m.cols());
88  bool status = invert(m,inv_m);
89  if (!status) {
90 #ifndef TO_CLEAN
91  std::ofstream out ("inv_failed.mtx");
93  out.close();
94 #endif // TO_CLEAN
95  error_macro("singular element matrix founded");
96  }
97  m = inv_m;
98 }
99 // ----------------------------------------------------------------------------
100 // instanciation in library
101 // ----------------------------------------------------------------------------
102 #define _RHEOLEF_instanciate(T) \
103 template T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m); \
104 template bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, const T&); \
105 template void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&); \
106 template void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, bool); \
107 
109 
110 }} // namespace rheolef::details
see the Float page for the full documentation
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
bool check_is_symmetric(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, const T &tol_m_max)
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
This file is part of Rheolef.
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:78
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
Float epsilon