59 DMType dm_name =
"DMMOFEM";
61 DMType dm_name_mg =
"DMMOFEM_MG";
69 auto core_log = logging::core::get();
71 LogManager::createSink(LogManager::getStrmWorld(),
"AT"));
73 LogManager::createSink(LogManager::getStrmWorld(),
"TIMER"));
74 LogManager::setLog(
"AT");
75 LogManager::setLog(
"TIMER");
85 simple->getAddBoundaryFE() =
true;
88 auto add_shared_entities_on_skeleton = [&]() {
90 auto boundary_meshset =
simple->getBoundaryMeshSet();
91 auto skeleton_meshset =
simple->getSkeletonMeshSet();
97 0,
simple->getDim() - 1, skeleton_ents,
true);
98 skeleton_ents = subtract(skeleton_ents, bdy_ents);
100 CHKERR m_field.
get_moab().add_entities(skeleton_meshset, skeleton_ents);
104 CHKERR add_shared_entities_on_skeleton();
114 const char *list_bases[] = {
"ainsworth",
"ainsworth_lobatto",
"demkowicz",
117 PetscInt choice_base_value = AINSWORTH;
119 LASBASETOP, &choice_base_value, &flg);
121 if (flg != PETSC_TRUE)
124 if (choice_base_value == AINSWORTH)
126 if (choice_base_value == AINSWORTH_LOBATTO)
128 else if (choice_base_value == DEMKOWICZ)
130 else if (choice_base_value == BERNSTEIN)
133 enum spaces { hdiv, hcurl, last_space };
134 const char *list_spaces[] = {
"hdiv",
"hcurl"};
135 PetscInt choice_space_value = hdiv;
137 last_space, &choice_space_value, &flg);
138 if (flg != PETSC_TRUE)
141 if (choice_space_value == hdiv)
143 else if (choice_space_value == hcurl)
149 CHKERR simple->addDomainBrokenField(
"BROKEN", space, base, 1);
160 CHKERR bc_mng->removeSideDOFs(
simple->getProblemName(),
"ZERO_FLUX",
165 auto assemble_domain_lhs = [&](
auto &pip) {
172 IT>::OpMixDivTimesScalar<SPACE_DIM>;
178 pip.push_back(
new OpHdivHdiv(
"BROKEN",
"BROKEN", beta));
179 auto unity = []() constexpr {
return 1; };
180 pip.push_back(
new OpHdivU(
"BROKEN",
"U", unity,
true));
188 op_loop_skeleton_side->getOpPtrVector(), {});
192 auto broken_data_ptr =
193 boost::make_shared<std::vector<BrokenBaseSideData>>();
198 op_loop_domain_side->getOpPtrVector(), {HDIV});
199 op_loop_domain_side->getOpPtrVector().push_back(
202 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
204 IT>::OpBrokenSpaceConstrain<1>;
205 op_loop_skeleton_side->getOpPtrVector().push_back(
206 new OpC(
"HYBRID", broken_data_ptr, 1.,
true,
false));
210 constexpr
int partition = 1;
212 op_print->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
216 if (
auto op_ptr =
dynamic_cast<BdyEleOp *
>(base_op_ptr)) {
217 auto fe_method = op_ptr->getFEMethod();
218 auto num_fe = fe_method->numeredEntFiniteElementPtr;
221 if (num_fe->getPStatus() & PSTATUS_SHARED)
222 MOFEM_LOG(
"SELF", Sev::inform) <<
"Num FE: " << *num_fe;
227 op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
230 pip.push_back(op_loop_skeleton_side);
235 auto assemble_domain_rhs = [&](
auto &pip) {
239 AT>::LinearForm<IT>::OpSource<1, 1>;
240 auto source = [&](
const double x,
const double y,
241 const double z) constexpr {
250 CHKERR assemble_domain_lhs(pip_mng->getOpDomainLhsPipeline());
251 CHKERR assemble_domain_rhs(pip_mng->getOpDomainRhsPipeline());
258 TetPolynomialBase::switchCacheBaseOn<HDIV>(
259 {pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
260 TetPolynomialBase::switchCacheBaseOn<L2>(
261 {pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
267 auto ksp = pip_mng->createKSP();
269 CHKERR KSPSetFromOptions(ksp);
270 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
271 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp";
273 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp <= Done";
275 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve";
277 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve <= Done";
279 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
280 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
284 auto ksp = pip_mng->createKSP();
286 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
287 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp";
288 CHKERR schur_ptr->setUp(ksp);
289 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp <= Done";
291 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve";
293 MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve <= Done";
295 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
296 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
301 auto check_residual = [&](
auto x,
auto f) {
307 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
313 auto div_flux_ptr = boost::make_shared<VectorDouble>();
315 "BROKEN", div_flux_ptr));
317 AT>::LinearForm<IT>::OpBaseTimesScalarField<1>;
319 domain_rhs.push_back(
new OpUDivFlux(
"U", div_flux_ptr, beta));
320 auto source = [&](
const double x,
const double y,
321 const double z) constexpr {
return 1; };
323 AT>::LinearForm<IT>::OpSource<1, 1>;
327 IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
329 AT>::LinearForm<IT>::OpBaseTimesVector<3, 3, 1>;
330 auto flux_ptr = boost::make_shared<MatrixDouble>();
331 domain_rhs.push_back(
333 boost::shared_ptr<VectorDouble> u_ptr =
334 boost::make_shared<VectorDouble>();
336 domain_rhs.push_back(
new OpHDivH(
"BROKEN", u_ptr, beta));
337 domain_rhs.push_back(
new OpHdivFlux(
"BROKEN", flux_ptr, beta));
345 op_loop_skeleton_side->getOpPtrVector(), {});
349 auto broken_data_ptr =
350 boost::make_shared<std::vector<BrokenBaseSideData>>();
355 op_loop_domain_side->getOpPtrVector(), {HDIV});
356 op_loop_domain_side->getOpPtrVector().push_back(
358 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
359 op_loop_domain_side->getOpPtrVector().push_back(
361 op_loop_domain_side->getOpPtrVector().push_back(
365 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
367 IT>::OpBrokenSpaceConstrainDHybrid<1>;
369 IT>::OpBrokenSpaceConstrainDFlux<1>;
370 op_loop_skeleton_side->getOpPtrVector().push_back(
371 new OpC_dHybrid(
"HYBRID", broken_data_ptr, 1.));
372 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
373 op_loop_skeleton_side->getOpPtrVector().push_back(
375 op_loop_skeleton_side->getOpPtrVector().push_back(
376 new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
379 domain_rhs.push_back(op_loop_skeleton_side);
382 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
383 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
385 pip_mng->getDomainRhsFE()->f =
f;
386 pip_mng->getSkeletonRhsFE()->f =
f;
387 pip_mng->getDomainRhsFE()->x = x;
388 pip_mng->getSkeletonRhsFE()->x = x;
391 simple->getDomainFEName(),
392 pip_mng->getDomainRhsFE());
394 CHKERR VecGhostUpdateBegin(
f, ADD_VALUES, SCATTER_REVERSE);
395 CHKERR VecGhostUpdateEnd(
f, ADD_VALUES, SCATTER_REVERSE);
400 CHKERR VecNorm(
f, NORM_2, &fnrm);
401 MOFEM_LOG_C(
"AT", Sev::inform,
"Residual %3.4e", fnrm);
403 constexpr
double eps = 1e-8;
406 "Residual norm larger than accepted");
411 auto calculate_error = [&]() {
415 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
421 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
422 auto flux_val_ptr = boost::make_shared<MatrixDouble>();
423 auto div_val_ptr = boost::make_shared<VectorDouble>();
424 auto source_ptr = boost::make_shared<VectorDouble>();
426 domain_rhs.push_back(
428 domain_rhs.push_back(
431 "BROKEN", div_val_ptr));
432 auto source = [&](
const double x,
const double y,
433 const double z) constexpr {
return -1; };
436 enum { DIV, GRAD,
LAST };
439 domain_rhs.push_back(
442 u_grad_ptr, mpi_vec, GRAD, flux_val_ptr));
445 simple->getDomainFEName(),
446 pip_mng->getDomainRhsFE());
447 CHKERR VecAssemblyBegin(mpi_vec);
448 CHKERR VecAssemblyEnd(mpi_vec);
451 const double *error_ind;
452 CHKERR VecGetArrayRead(mpi_vec, &error_ind);
454 <<
"Approximation error ||div flux - source||: "
455 << std::sqrt(error_ind[DIV]);
456 MOFEM_LOG(
"AT", Sev::inform) <<
"Approximation error ||grad-flux||: "
457 << std::sqrt(error_ind[GRAD]);
458 CHKERR VecRestoreArrayRead(mpi_vec, &error_ind);
464 auto get_post_proc_fe = [&]() {
467 auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
473 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
476 op_loop_side->getOpPtrVector(), {HDIV});
477 auto u_vec_ptr = boost::make_shared<VectorDouble>();
478 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
479 op_loop_side->getOpPtrVector().push_back(
481 op_loop_side->getOpPtrVector().push_back(
483 op_loop_side->getOpPtrVector().push_back(
487 post_proc_fe->getPostProcMesh(),
489 post_proc_fe->getMapGaussPts(),
493 {{
"BROKEN", flux_mat_ptr}},
502 auto post_proc_fe = get_post_proc_fe();
504 simple->getBoundaryFEName(), post_proc_fe);
505 CHKERR post_proc_fe->writeFile(
"out_result.h5m");