11 static char help[] =
"...\n\n";
30 constexpr
double a = 1;
31 constexpr
double a2 =
a *
a;
34 constexpr
double A = 6371220;
40 auto res_J = [](
const double x,
const double y,
const double z) {
41 const double res = (x * x + y * y + z * z -
a2);
45 auto res_J_dx = [](
const double x,
const double y,
const double z) {
46 const double res =
res_J(x, y, z);
51 auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
52 const double res =
res_J(x, y, z);
55 (res * 2 + (4 * x * x)),
60 (2 * res + (4 * y * y)),
65 (2 * res + (4 * z * z))};
71 boost::shared_ptr<MatrixDouble> dot_x_ptr)
73 xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
78 auto t_w = getFTensor0IntegrationWeight();
81 auto t_x0 = getFTensor1CoordsAtGaussPts();
82 auto t_x = getFTensor1FromMat<3>(*xPtr);
83 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
84 auto t_normal = getFTensor1NormalsAtGaussPts();
88 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
94 t_P(
i,
j) = t_n(
i) * t_n(
j);
97 auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
99 const double alpha = t_w;
100 auto t_nf = getFTensor1FromArray<3, 3>(locF);
101 double l = std::sqrt(t_normal(
i) * t_normal(
i));
105 alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
108 for (; rr != nbRows / 3; ++rr) {
110 t_nf(
j) += t_row_base * t_res(
j);
115 for (; rr < nbRowBaseFunctions; ++rr) {
130 boost::shared_ptr<MatrixDouble>
xPtr;
137 boost::shared_ptr<MatrixDouble> dot_x_ptr)
140 xPtr(x_ptr), xDotPtr(dot_x_ptr) {
148 auto t_w = getFTensor0IntegrationWeight();
151 auto t_x0 = getFTensor1CoordsAtGaussPts();
152 auto t_x = getFTensor1FromMat<3>(*xPtr);
153 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
154 auto t_normal = getFTensor1NormalsAtGaussPts();
156 auto get_t_mat = [&](
const int rr) {
158 &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
160 &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
162 &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
166 const double ts_a = getTSa();
168 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
174 t_P(
i,
j) = t_n(
i) * t_n(
j);
177 auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
178 double l = std::sqrt(t_normal(
i) * t_normal(
i));
180 const double alpha = t_w;
183 (alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
186 for (; rr != nbRows / 3; rr++) {
189 auto t_mat = get_t_mat(3 * rr);
191 for (
int cc = 0; cc != nbCols / 3; cc++) {
193 const double rc = t_row_base * t_col_base;
194 t_mat(
i,
j) += rc * t_lhs(
i,
j);
202 for (; rr < nbRowBaseFunctions; ++rr)
215 boost::shared_ptr<MatrixDouble>
xPtr;
225 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
232 auto t_w = getFTensor0IntegrationWeight();
233 auto t_x = getFTensor1FromMat<3>(*xPtr);
234 auto t_normal = getFTensor1NormalsAtGaussPts();
235 auto nb_integration_pts = getGaussPts().size2();
239 for (
int gg = 0; gg != nb_integration_pts; gg++) {
241 double l = std::sqrt(t_normal(
i) * t_normal(
i));
242 error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
249 CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
257 boost::shared_ptr<MatrixDouble>
xPtr;
335 auto x_ptr = boost::make_shared<MatrixDouble>();
336 auto dot_x_ptr = boost::make_shared<MatrixDouble>();
337 auto det_ptr = boost::make_shared<VectorDouble>();
338 auto jac_ptr = boost::make_shared<MatrixDouble>();
339 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
341 auto def_ops = [&](
auto &pipeline) {
352 new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
354 new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
366 CHKERR mField.loop_dofs(
"HO_POSITIONS", ent_method_material);
371 auto dm =
simple->getDM();
373 ts = pipeline_mng->createTSIM();
376 CHKERR TSSetMaxSteps(ts, 1);
377 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
382 CHKERR TSSetSolution(ts, T);
383 CHKERR TSSetFromOptions(ts);
386 CHKERR TSGetTime(ts, &ftime);
397 auto x_ptr = boost::make_shared<MatrixDouble>();
398 auto det_ptr = boost::make_shared<VectorDouble>();
399 auto jac_ptr = boost::make_shared<MatrixDouble>();
400 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
403 auto dm =
simple->getDM();
406 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
408 post_proc_fe->getOpPtrVector().push_back(
413 post_proc_fe->getOpPtrVector().push_back(
415 new OpPPMap(post_proc_fe->getPostProcMesh(),
416 post_proc_fe->getMapGaussPts(),
420 {{
"HO_POSITIONS", x_ptr}},
429 CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
431 auto error_fe = boost::make_shared<DomainEle>(mField);
433 error_fe->getOpPtrVector().push_back(
435 error_fe->getOpPtrVector().push_back(
437 error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
439 error_fe->preProcessHook = [&]() {
441 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
443 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
448 error_fe->postProcessHook = [&]() {
455 <<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
464 if (mField.get_comm_size() > 1)
465 CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m",
"MOAB",
466 "PARALLEL=WRITE_PART");
468 CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m");
473 int main(
int argc,
char *argv[]) {
476 const char param_file[] =
"param_file.petsc";
479 auto core_log = logging::core::get();
488 DMType dm_name =
"DMMOFEM";