Rheolef  7.1
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 {
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  check_macro (invert(m,inv_m), "singular element matrix founded");
89  m = inv_m;
90 }
91 // ----------------------------------------------------------------------------
92 // instanciation in library
93 // ----------------------------------------------------------------------------
94 #define _RHEOLEF_instanciate(T) \
95 template T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m); \
96 template bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, const T&); \
97 template void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&); \
98 template void local_invert (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, bool); \
99 
101 
102 }// namespace rheolef
rheolef::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::local_invert
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
Definition: form_vf_assembly.cc:77
rheolef::local_lump
void local_lump(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: form_vf_assembly.cc:61
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::norm_max
T norm_max(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: form_vf_assembly.cc:33
Float
see the Float page for the full documentation
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
rheolef::_RHEOLEF_instanciate
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
epsilon
Float epsilon
Definition: transmission_error.cc:25
T
Expr1::float_type T
Definition: field_expr.h:261