static char help[] =
"testing field evaluator\n\n";
template <
typename OP>
struct MyOp :
public OP {
MatrixShallowArrayAdaptor<double> evalPoints;
MyOp(std::array<double, 12> &eval_points)
:
OP(
"FIELD1", OP::OPROW),
evalPoints(
MOFEM_LOG(
"SELF", Sev::inform) <<
"FE " << OP::getFEEntityHandle();
MOFEM_LOG(
"SELF", Sev::inform) <<
"Integration pts" << std::endl;
MOFEM_LOG(
"SELF", Sev::inform) << OP::getGaussPts() << endl;
MOFEM_LOG(
"SELF", Sev::inform) <<
"Global coordinates " << endl;
MOFEM_LOG(
"SELF", Sev::inform) << OP::getCoordsAtGaussPts() << std::endl;
for (int gg = 0; gg != OP::getCoordsAtGaussPts().size1(); ++gg) {
int pt_number = OP::getGaussPts()(OP::getGaussPts().size1() - 1, gg);
MOFEM_LOG(
"SELF", Sev::inform) <<
"gg " << gg << std::endl;
MOFEM_LOG(
"SELF", Sev::inform) <<
"pt " << pt_number << std::endl;
ublas::matrix_row<MatrixDouble> coord_at_gauss_pt(
OP::getCoordsAtGaussPts(), gg);
ublas::matrix_row<MatrixShallowArrayAdaptor<double>> eval_coord(
evalPoints, pt_number);
MOFEM_LOG(
"SELF", Sev::inform) <<
"coord_at_gauss_pt ";
MOFEM_LOG(
"SELF", Sev::inform) << coord_at_gauss_pt << std::endl;
MOFEM_LOG(
"SELF", Sev::inform) <<
"eval_coord ";
MOFEM_LOG(
"SELF", Sev::inform) << eval_coord << std::endl;
double error = norm_2(coord_at_gauss_pt - eval_coord);
if (error > 1e-12)
"Difference at %d error = %3.4e", pt_number, error);
}
}
}
};
auto get_rule = [](
int order_row,
int order_col,
int order_data) {
return -1;
};
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
{
auto dm = simple_interface->
getDM();
if (simple_interface->
getDim() == 3) {
std::array<double, 3> point = {0, 0, 0};
std::array<double, 12> eval_points = {0.0, 0.0, 0.0, 0.0, 0.01, 0.0,
-0.1, 0.0, 0.0, 0.1, 0.0, 0.0};
if (auto fe_method = data->feMethodPtr.lock()) {
fe_method->getOpPtrVector().push_back(
new MyOp<VolOp>(eval_points));
} else
"Pointer to element does not exists");
data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
}
if (simple_interface->
getDim() == 2) {
std::array<double, 3> point = {0, 0, 0};
std::array<double, 12> eval_points = {
0.0, 0.0, 0.0,
0.0, 0.01, 0.0,
-0.1, 0.0, 0.0,
0.1, 0.0, 0.0};
if (auto fe_method = data->feMethodPtr.lock()) {
fe_method->getOpPtrVector().push_back(
new MyOp<FaceOp>(eval_points));
} else
"Pointer to element does not exists");
data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
}
}
}
return 0;
}