Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
form_vf_assembly.cc
Go to the documentation of this file.
1
21//
22#include "rheolef/form_vf_assembly.h"
23#include "rheolef/eigen_util.h"
24
25namespace rheolef { namespace details {
26using namespace std;
27
28// ------------------------------------------------------------------------------
29// norm_max
30// ------------------------------------------------------------------------------
31template <class T>
32T
33norm_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// ------------------------------------------------------------------------------
44template <class T>
45bool
46check_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// ------------------------------------------------------------------------------
59template <class T>
60void
61local_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// ------------------------------------------------------------------------------
75template <class T>
76void
77local_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");
92 put_matrix_market (out, m);
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) \
103template T norm_max (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m); \
104template bool check_is_symmetric (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&, const T&); \
105template void local_lump (Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&); \
106template 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 _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
#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.
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition eigen_util.h:78
STL namespace.