v0.9.2
PoissonOperators.hpp
#ifndef __POISSONOPERTATORS_HPP__
#define __POISSONOPERTATORS_HPP__
#include <stdlib.h>
namespace PoissonOps {
// Operators
const double body_source = 1;
struct OpLhsUU : public OpVolEle {
public:
OpLhsUU(std::string domain_field1, std::string domain_field2)
: OpVolEle(domain_field1, domain_field2, OpVolEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locLhs.resize(nb_row_dofs, nb_col_dofs, false);
locLhs.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = t_w * vol;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; ++cc) {
locLhs(rr, cc) += a * t_row_diff_base(i) * t_col_diff_base(i);
++t_col_diff_base;
} // endFor cc
++t_row_diff_base;
} // endFor rr
++t_w;
} // endFor gg
CHKERR MatSetValues(getFEMethod()->ksp_B, row_data, col_data,
&locLhs(0, 0), ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues(getFEMethod()->ksp_B, col_data, row_data,
&transLocLhs(0, 0), ADD_VALUES);
} // (row_side != col_side || row_type != col_type)
} // endIf (nb_row_dofs && nb_col_dofs)
}
private:
};
struct OpRhsU : public OpVolEle {
public:
OpRhsU(std::string domain_field) : OpVolEle(domain_field, OpVolEle::OPROW) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
locRhs.resize(nb_dofs, false);
locRhs.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_base = data.getFTensor0N();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
for (int rr = 0; rr != nb_dofs; ++rr) {
locRhs[rr] += a * (t_base * body_source);
++t_base;
} // endFor rr
++t_w;
} // endFor gg
CHKERR VecSetOption(getFEMethod()->ksp_f, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
CHKERR VecSetValues(getFEMethod()->ksp_f, data, &*locRhs.begin(),
ADD_VALUES);
} // endIf (nb_dofs)
}
private:
};
}; // namespace PoissonOps
#endif //__POISSONOPERTATORS_HPP__