1 #ifndef _RHEO_FORM_ELEMENT_H
2 #define _RHEO_FORM_ELEMENT_H
31 #include "rheolef/quadrature.h"
32 #include "rheolef/basis.h"
46 template<
class T,
class Value,
class Basis>
51 Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
59 Eigen::Matrix<Value,Eigen::Dynamic,1>
phi (loc_ndof);
60 mass.resize (loc_ndof, loc_ndof);
62 last_q = quad.
end(hat_K); iter_q != last_q; iter_q++) {
63 b.evaluate (hat_K, (*iter_q).x,
phi);
64 for (
size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
65 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
67 if (fabs(coeff_w) < eps)
continue;
68 mass.coeffRef(loc_idof,loc_jdof) += coeff_w;
73 mass.makeCompressed();
75 template<
class T,
class Basis>
80 Eigen::SparseMatrix<T,Eigen::RowMajor>& mass)
82 switch (
b.valued_tag()) {
85 default:
error_macro(
"tensor-like mass_element test: not yet");
91 template<
class T,
class Basis>
96 Eigen::SparseMatrix<T,Eigen::RowMajor>& mass_bdr)
104 Eigen::Matrix<T,Eigen::Dynamic,1>
phi (loc_ndof);
105 mass_bdr.resize (loc_ndof, loc_ndof);
106 for (
size_type loc_isid = 0, loc_nsid = tilde_K.
n_side(); loc_isid < loc_nsid; ++loc_isid) {
111 last_q = quad.
end(hat_S); iter_q != last_q; iter_q++) {
112 T tilde_wq = tilde_Js*(*iter_q).w;
113 if (
b.option().is_restricted_to_sides()) {
114 b.evaluate_on_side (tilde_K, sid, (*iter_q).x,
phi);
116 b.evaluate (tilde_K, (*iter_q).x,
phi);
118 for (
size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
119 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
121 if (fabs(coeff_w) < eps)
continue;
122 mass_bdr.coeffRef(loc_idof,loc_jdof) += coeff_w;
130 template<
class T,
class Basis>
135 Eigen::SparseMatrix<T,Eigen::RowMajor>& grad_grad)
143 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi (loc_ndof);
144 grad_grad.resize (loc_ndof, loc_ndof);
146 last_q = quad.
end(hat_K); iter_q != last_q; iter_q++) {
147 b.grad_evaluate (hat_K, (*iter_q).x, grad_phi);
148 for (
size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
149 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
150 T coeff_w = (*iter_q).w*
dot(grad_phi[loc_idof],grad_phi[loc_jdof]);
151 if (fabs(coeff_w) < eps)
continue;
152 grad_grad.coeffRef (loc_idof,loc_jdof) += coeff_w;
158 #endif // _RHEO_FORM_ELEMENT_H