static char help[] =
"...\n\n";
};
constexpr
double a2 =
a *
a;
constexpr
double a4 =
a2 *
a2;
constexpr
double A = 6371220;
auto res_J = [](
const double x,
const double y,
const double z) {
const double res = (x * x + y * y + z * z -
a2);
return res;
};
auto res_J_dx = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
res * (2 * z)};
};
auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
(res * 2 + (4 * x * x)),
(4 * y * x),
(4 * z * x),
(4 * x * y),
(2 * res + (4 * y * y)),
(4 * z * y),
(4 * x * z),
(4 * y * z),
(2 * res + (4 * z * z))};
};
boost::shared_ptr<MatrixDouble> dot_x_ptr)
xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
auto t_w = getFTensor0IntegrationWeight();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
const double alpha = t_w;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
double l = std::sqrt(t_normal(
i) * t_normal(
i));
alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(
j) += t_row_base * t_res(
j);
++t_row_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
}
++t_w;
++t_x;
++t_dot_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
boost::shared_ptr<MatrixDouble> dot_x_ptr)
xPtr(x_ptr), xDotPtr(dot_x_ptr) {
this->sYmm = false;
}
auto t_w = getFTensor0IntegrationWeight();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const double ts_a = getTSa();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
double l = std::sqrt(t_normal(
i) * t_normal(
i));
const double alpha = t_w;
(alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
const double rc = t_row_base * t_col_base;
t_mat(
i,
j) += rc * t_lhs(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
xPtr(x_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
auto t_w = getFTensor0IntegrationWeight();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto nb_integration_pts = getGaussPts().size2();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(
i) * t_normal(
i));
error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
++t_w;
++t_x;
++t_normal;
}
CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
private:
};
}
}
}
}
};
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
pipeline.push_back(
};
new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
}
CHKERR mField.loop_dofs(
"HO_POSITIONS", ent_method_material);
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
}
auto x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{},
{{"HO_POSITIONS", x_ptr}},
{}, {}
)
);
CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(mField);
error_fe->getOpPtrVector().push_back(
error_fe->getOpPtrVector().push_back(
error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
error_fe->preProcessHook = [&]() {
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
};
error_fe->postProcessHook = [&]() {
double error2;
<<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
};
if (mField.get_comm_size() > 1)
CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m",
"MOAB",
"PARALLEL=WRITE_PART");
else
CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m");
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}