v0.15.0
Loading...
Searching...
No Matches
test_broken_space.cpp
Go to the documentation of this file.
1/**
2 * \example test_broken_space.cpp
3 *
4 * Testing broken spaces. Implementations works for 2d and 3d meshes, is aiming
5 * to test H-div broken base functions, and L2 base on skeleton.
6 *
7 * Also, it test block matrix with Schur complement.
8 *
9 */
10
11#include <MoFEM.hpp>
13
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
18constexpr bool debug = false;
19
20constexpr AssemblyType BlockAssemblyType = (BLOCK_ASSEMBLE_SELECTION)
21 ? AssemblyType::BLOCK_MAT
23
24constexpr AssemblyType AT =
26 : AssemblyType::PETSC; //< selected assembly type
27
29 IntegrationType::GAUSS; //< selected integration type
30
31constexpr int SPACE_DIM =
32 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
33
38
40using DomainEleOp = DomainEle::UserDataOperator;
41using BdyEleOp = BoundaryEle::UserDataOperator;
42using SideEleOp = EleOnSide::UserDataOperator;
43
44struct SetUpSchur {
45 static boost::shared_ptr<SetUpSchur>
48
49protected:
50 SetUpSchur() = default;
51 virtual ~SetUpSchur() = default;
52};
53
55
56int main(int argc, char *argv[]) {
57
58 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
59
60 try {
61
62 //! [Register MoFEM discrete manager in PETSc]
63 DMType dm_name = "DMMOFEM";
64 CHKERR DMRegister_MoFEM(dm_name);
65 DMType dm_name_mg = "DMMOFEM_MG";
67 //! [Register MoFEM discrete manager in PETSc
68
69 moab::Core mb_instance;
70 moab::Interface &moab = mb_instance;
71
72 // Add logging channel for example
73 auto core_log = logging::core::get();
74 core_log->add_sink(
76 core_log->add_sink(
79 LogManager::setLog("TIMER");
80 MOFEM_LOG_TAG("AT", "atom_test");
81 MOFEM_LOG_TAG("TIMER", "timer");
82
83 // Create MoFEM instance
84 MoFEM::Core core(moab);
85 MoFEM::Interface &m_field = core;
86
87 auto *simple = m_field.getInterface<Simple>();
89 simple->getAddBoundaryFE() = true;
90 CHKERR simple->loadFile();
91
92 auto add_shared_entities_on_skeleton = [&]() {
94 auto boundary_meshset = simple->getBoundaryMeshSet();
95 auto skeleton_meshset = simple->getSkeletonMeshSet();
96 Range bdy_ents;
97 CHKERR m_field.get_moab().get_entities_by_handle(boundary_meshset,
98 bdy_ents, true);
99 Range skeleton_ents;
100 CHKERR m_field.get_moab().get_entities_by_dimension(
101 0, simple->getDim() - 1, skeleton_ents, true);
102 skeleton_ents = subtract(skeleton_ents, bdy_ents);
103 CHKERR m_field.get_moab().clear_meshset(&skeleton_meshset, 1);
104 CHKERR m_field.get_moab().add_entities(skeleton_meshset, skeleton_ents);
106 };
107
108 CHKERR add_shared_entities_on_skeleton();
109
110 // Declare elements
111 enum bases {
112 AINSWORTH,
113 AINSWORTH_LOBATTO,
114 DEMKOWICZ,
115 BERNSTEIN,
116 LASBASETOP
117 };
118 const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
119 "bernstein"};
120 PetscBool flg;
121 PetscInt choice_base_value = AINSWORTH;
122 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
123 LASBASETOP, &choice_base_value, &flg);
124
125 if (flg != PETSC_TRUE)
126 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
128 if (choice_base_value == AINSWORTH)
130 if (choice_base_value == AINSWORTH_LOBATTO)
132 else if (choice_base_value == DEMKOWICZ)
134 else if (choice_base_value == BERNSTEIN)
136
137 enum spaces { hdiv, hcurl, last_space };
138 const char *list_spaces[] = {"hdiv", "hcurl"};
139 PetscInt choice_space_value = hdiv;
140 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-space", list_spaces,
141 last_space, &choice_space_value, &flg);
142 if (flg != PETSC_TRUE)
143 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
144 FieldSpace space = HDIV;
145 if (choice_space_value == hdiv)
146 space = HDIV;
147 else if (choice_space_value == hcurl)
148 space = HCURL;
149
150 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &approx_order,
151 PETSC_NULLPTR);
152
153 CHKERR simple->addDomainBrokenField("BROKEN", space, base, 1);
154 CHKERR simple->addDomainField("U", L2, base, 1);
155 CHKERR simple->addSkeletonField("HYBRID", L2, base, 1);
156
157 CHKERR simple->setFieldOrder("BROKEN", approx_order);
158 CHKERR simple->setFieldOrder("U", approx_order - 1);
159 CHKERR simple->setFieldOrder("HYBRID", approx_order - 1);
160
161 CHKERR simple->setUp();
162
163 auto bc_mng = m_field.getInterface<BcManager>();
164 CHKERR bc_mng->removeSideDOFs(simple->getProblemName(), "ZERO_FLUX",
165 "BROKEN", SPACE_DIM, 0, 1, true);
166
167 auto integration_rule = [](int, int, int p) { return 2 * p; };
168
169 auto assemble_domain_lhs = [&](auto &pip) {
172
176 IT>::OpMixDivTimesScalar<SPACE_DIM>;
177
178 auto beta = [](const double, const double, const double) constexpr {
179 return 1;
180 };
181
182 pip.push_back(new OpHdivHdiv("BROKEN", "BROKEN", beta));
183 auto unity = []() constexpr { return 1; };
184 pip.push_back(new OpHdivU("BROKEN", "U", unity, true));
185
186 // First: Iterate over skeleton FEs adjacent to Domain FEs
187 // Note: BoundaryEle, i.e. uses skeleton interation rule
188 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
189 m_field, simple->getSkeletonFEName(), SPACE_DIM - 1, Sev::noisy);
190 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = integration_rule;
192 op_loop_skeleton_side->getOpPtrVector(), {});
193
194 // Second: Iterate over domain FEs adjacent to skelton, particularly one
195 // domain element.
196 auto broken_data_ptr =
197 boost::make_shared<std::vector<BrokenBaseSideData>>();
198 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
199 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
200 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy);
202 op_loop_domain_side->getOpPtrVector(), {HDIV});
203 op_loop_domain_side->getOpPtrVector().push_back(
204 new OpGetBrokenBaseSideData<SideEleOp>("BROKEN", broken_data_ptr));
205
206 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
208 IT>::OpBrokenSpaceConstrain<1>;
209 op_loop_skeleton_side->getOpPtrVector().push_back(
210 new OpC("HYBRID", broken_data_ptr, 1., true, false));
211
212 if (debug) {
213 // print skeleton elements on partition
214 constexpr int partition = 1;
215 auto op_print = new BdyEleOp(NOSPACE, BdyEleOp::OPSPACE);
216 op_print->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
217 EntityType type,
220 if (auto op_ptr = dynamic_cast<BdyEleOp *>(base_op_ptr)) {
221 auto fe_method = op_ptr->getFEMethod();
222 auto num_fe = fe_method->numeredEntFiniteElementPtr;
223
224 if (m_field.get_comm_rank() == partition) {
225 if (num_fe->getPStatus() & PSTATUS_SHARED)
226 MOFEM_LOG("SELF", Sev::inform) << "Num FE: " << *num_fe;
227 }
228 }
230 };
231 op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
232 };
233
234 pip.push_back(op_loop_skeleton_side);
235
237 };
238
239 auto assemble_domain_rhs = [&](auto &pip) {
244 auto source = [&](const double x, const double y,
245 const double z) constexpr {
246 return -1; // sin(100 * (x / 10.) * M_PI_2);
247 };
248 pip.push_back(new OpDomainSource("U", source));
250 };
251
252 auto *pip_mng = m_field.getInterface<PipelineManager>();
253
254 CHKERR assemble_domain_lhs(pip_mng->getOpDomainLhsPipeline());
255 CHKERR assemble_domain_rhs(pip_mng->getOpDomainRhsPipeline());
256
257 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
258 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
259 CHKERR pip_mng->setSkeletonLhsIntegrationRule(integration_rule);
260 CHKERR pip_mng->setSkeletonRhsIntegrationRule(integration_rule);
261
262 TetPolynomialBase::switchCacheBaseOn<HDIV>(
263 {pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
264 TetPolynomialBase::switchCacheBaseOn<L2>(
265 {pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
266
267 auto x = createDMVector(simple->getDM());
268 auto f = vectorDuplicate(x);
269
270 if (AT == PETSC) {
271 auto ksp = pip_mng->createKSP();
272
273 CHKERR KSPSetFromOptions(ksp);
274 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
275 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
276 CHKERR KSPSetUp(ksp);
277 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
278
279 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
280 CHKERR KSPSolve(ksp, f, x);
281 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
282
283 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
284 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
285 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
286 SCATTER_REVERSE);
287 } else {
288 auto ksp = pip_mng->createKSP();
289 auto schur_ptr = SetUpSchur::createSetUpSchur(m_field);
290 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
291 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
292 CHKERR schur_ptr->setUp(ksp);
293 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
294
295 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
296 CHKERR KSPSolve(ksp, f, x);
297 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
298
299 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
300 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
301 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
302 SCATTER_REVERSE);
303 }
304
305 auto check_residual = [&](auto x, auto f) {
307 auto *simple = m_field.getInterface<Simple>();
308 auto *pip_mng = m_field.getInterface<PipelineManager>();
309
310 // auto &skeleton_rhs = pip_mng->getOpSkeletonRhsPipeline();
311 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
312 // skeleton_rhs.clear();
313 domain_rhs.clear();
314
316
317 auto div_flux_ptr = boost::make_shared<VectorDouble>();
318 domain_rhs.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
319 "BROKEN", div_flux_ptr));
322 auto beta = [](double, double, double) constexpr { return 1; };
323 domain_rhs.push_back(new OpUDivFlux("U", div_flux_ptr, beta));
324 auto source = [&](const double x, const double y,
325 const double z) constexpr { return 1; };
328 domain_rhs.push_back(new OpDomainSource("U", source));
329
331 IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
334 auto flux_ptr = boost::make_shared<MatrixDouble>();
335 domain_rhs.push_back(
336 new OpCalculateHVecVectorField<3>("BROKEN", flux_ptr));
337 boost::shared_ptr<VectorDouble> u_ptr =
338 boost::make_shared<VectorDouble>();
339 domain_rhs.push_back(new OpCalculateScalarFieldValues("U", u_ptr));
340 domain_rhs.push_back(new OpHDivH("BROKEN", u_ptr, beta));
341 domain_rhs.push_back(new OpHdivFlux("BROKEN", flux_ptr, beta));
342
343 // First: Iterate over skeleton FEs adjacent to Domain FEs
344 // Note: BoundaryEle, i.e. uses skeleton interation rule
345 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
346 m_field, simple->getSkeletonFEName(), SPACE_DIM - 1, Sev::noisy);
347 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = integration_rule;
349 op_loop_skeleton_side->getOpPtrVector(), {});
350
351 // Second: Iterate over domain FEs adjacent to skelton, particularly one
352 // domain element.
353 auto broken_data_ptr =
354 boost::make_shared<std::vector<BrokenBaseSideData>>();
355 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
356 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
357 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy);
359 op_loop_domain_side->getOpPtrVector(), {HDIV});
360 op_loop_domain_side->getOpPtrVector().push_back(
361 new OpGetBrokenBaseSideData<SideEleOp>("BROKEN", broken_data_ptr));
362 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
363 op_loop_domain_side->getOpPtrVector().push_back(
364 new OpCalculateHVecTensorField<1, 3>("BROKEN", flux_mat_ptr));
365 op_loop_domain_side->getOpPtrVector().push_back(
366 new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
367
368 // Assemble on skeleton
369 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
371 IT>::OpBrokenSpaceConstrainDHybrid<1>;
373 IT>::OpBrokenSpaceConstrainDFlux<1>;
374 op_loop_skeleton_side->getOpPtrVector().push_back(
375 new OpC_dHybrid("HYBRID", broken_data_ptr, 1.));
376 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
377 op_loop_skeleton_side->getOpPtrVector().push_back(
378 new OpCalculateVectorFieldValues<1>("HYBRID", hybrid_ptr));
379 op_loop_skeleton_side->getOpPtrVector().push_back(
380 new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
381
382 // Add skeleton to domain pipeline
383 domain_rhs.push_back(op_loop_skeleton_side);
384
385 CHKERR VecZeroEntries(f);
386 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
387 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
388
389 pip_mng->getDomainRhsFE()->f = f;
390 pip_mng->getSkeletonRhsFE()->f = f;
391 pip_mng->getDomainRhsFE()->x = x;
392 pip_mng->getSkeletonRhsFE()->x = x;
393
395 simple->getDomainFEName(),
396 pip_mng->getDomainRhsFE());
397
398 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
399 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
400 CHKERR VecAssemblyBegin(f);
401 CHKERR VecAssemblyEnd(f);
402
403 double fnrm;
404 CHKERR VecNorm(f, NORM_2, &fnrm);
405 MOFEM_LOG_C("AT", Sev::inform, "Residual %3.4e", fnrm);
406
407 constexpr double eps = 1e-8;
408 if (fnrm > eps)
409 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
410 "Residual norm larger than accepted");
411
413 };
414
415 auto calculate_error = [&]() {
417
418 // auto &skeleton_rhs = pip_mng->getOpSkeletonRhsPipeline();
419 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
420 // skeleton_rhs.clear();
421 domain_rhs.clear();
422
424
425 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
426 auto flux_val_ptr = boost::make_shared<MatrixDouble>();
427 auto div_val_ptr = boost::make_shared<VectorDouble>();
428 auto source_ptr = boost::make_shared<VectorDouble>();
429
430 domain_rhs.push_back(
431 new OpCalculateScalarFieldGradient<SPACE_DIM>("U", u_grad_ptr));
432 domain_rhs.push_back(
433 new OpCalculateHVecVectorField<3, SPACE_DIM>("BROKEN", flux_val_ptr));
434 domain_rhs.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
435 "BROKEN", div_val_ptr));
436 auto source = [&](const double x, const double y,
437 const double z) constexpr { return -1; };
438 domain_rhs.push_back(new OpGetTensor0fromFunc(source_ptr, source));
439
440 enum { DIV, GRAD, LAST };
441 auto mpi_vec = createVectorMPI(
442 m_field.get_comm(), (!m_field.get_comm_rank()) ? LAST : 0, LAST);
443 domain_rhs.push_back(
444 new OpCalcNormL2Tensor0(div_val_ptr, mpi_vec, DIV, source_ptr));
445 domain_rhs.push_back(new OpCalcNormL2Tensor1<SPACE_DIM>(
446 u_grad_ptr, mpi_vec, GRAD, flux_val_ptr));
447
449 simple->getDomainFEName(),
450 pip_mng->getDomainRhsFE());
451 CHKERR VecAssemblyBegin(mpi_vec);
452 CHKERR VecAssemblyEnd(mpi_vec);
453
454 if (!m_field.get_comm_rank()) {
455 const double *error_ind;
456 CHKERR VecGetArrayRead(mpi_vec, &error_ind);
457 MOFEM_LOG("AT", Sev::inform)
458 << "Approximation error ||div flux - source||: "
459 << std::sqrt(error_ind[DIV]);
460 MOFEM_LOG("AT", Sev::inform) << "Approximation error ||grad-flux||: "
461 << std::sqrt(error_ind[GRAD]);
462 CHKERR VecRestoreArrayRead(mpi_vec, &error_ind);
463 }
464
466 };
467
468 auto get_post_proc_fe = [&]() {
471 auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
472
473 auto op_loop_side = new OpLoopSide<EleOnSide>(
474 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy,
475 boost::make_shared<
476 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
477 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
478
480 op_loop_side->getOpPtrVector(), {HDIV});
481 auto u_vec_ptr = boost::make_shared<VectorDouble>();
482 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
483 op_loop_side->getOpPtrVector().push_back(
484 new OpCalculateScalarFieldValues("U", u_vec_ptr));
485 op_loop_side->getOpPtrVector().push_back(
486 new OpCalculateHVecVectorField<3>("BROKEN", flux_mat_ptr));
487 op_loop_side->getOpPtrVector().push_back(
488
489 new OpPPMap(
490
491 post_proc_fe->getPostProcMesh(),
492
493 post_proc_fe->getMapGaussPts(),
494
495 {{"U", u_vec_ptr}},
496
497 {{"BROKEN", flux_mat_ptr}},
498
499 {}, {})
500
501 );
502
503 return post_proc_fe;
504 };
505
506 auto post_proc_fe = get_post_proc_fe();
508 simple->getBoundaryFEName(), post_proc_fe);
509 CHKERR post_proc_fe->writeFile("out_result.h5m");
510
511 CHKERR calculate_error();
512 CHKERR check_residual(x, f);
513 }
515
517}
518
519struct SetUpSchurImpl : public SetUpSchur {
520
522
523 virtual ~SetUpSchurImpl() = default;
524
526
527private:
530 };
531
535 auto pip_mng = mField.getInterface<PipelineManager>();
536
537 CHKERR KSPSetFromOptions(ksp);
538 PC pc;
539 CHKERR KSPGetPC(ksp, &pc);
540
541 PetscBool is_pcfs = PETSC_FALSE;
542 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
543 if (is_pcfs) {
544
545 MOFEM_LOG("AT", Sev::inform) << "Setup Schur pc";
546
547 auto create_sub_dm = [&]() {
549
550 auto create_dm = [&](
551
552 std::string problem_name,
553 std::vector<std::string> fe_names,
554 std::vector<std::string> fields,
555
556 auto dm_type
557
558 ) {
559 auto dm = createDM(mField.get_comm(), dm_type);
560 auto create_dm_imp = [&]() {
562 CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), problem_name.c_str());
563 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
564 for (auto fe : fe_names) {
566 }
567 CHKERR DMMoFEMAddElement(dm, simple->getSkeletonFEName());
568 for (auto field : fields) {
569 CHKERR DMMoFEMAddSubFieldRow(dm, field);
570 CHKERR DMMoFEMAddSubFieldCol(dm, field);
571 }
572 CHKERR DMSetUp(dm);
574 };
576 create_dm_imp(),
577 "Error in creating schurDM. It is possible that schurDM is "
578 "already created");
579 return dm;
580 };
581
582 auto schur_dm = create_dm(
583
584 "SCHUR",
585
586 {simple->getDomainFEName(), simple->getSkeletonFEName()},
587
588 {"HYBRID"},
589
590 "DMMOFEM_MG");
591
592 auto block_dm = create_dm(
593
594 "BLOCK",
595
596 {simple->getDomainFEName(), simple->getSkeletonFEName()},
597
598 {"BROKEN", "U"},
599
600 "DMMOFEM");
601
602 return std::make_tuple(schur_dm, block_dm);
603 };
604
605 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
606 auto block_mat_data = createBlockMatStructure(
607 simple->getDM(),
608
609 {
610
611 {
612
613 simple->getDomainFEName(),
614
615 {
616
617 {"BROKEN", "BROKEN"},
618 {"U", "U"},
619 {"BROKEN", "U"},
620 {"U", "BROKEN"}
621
622 }
623
624 },
625
626 {
627
628 simple->getSkeletonFEName(),
629
630 {
631
632 {"BROKEN", "HYBRID"}, {"HYBRID", "BROKEN"}
633
634 }
635
636 }
637
638 }
639
640 );
641
643
644 {schur_dm, block_dm}, block_mat_data,
645
646 {"BROKEN", "U"}, {nullptr, nullptr}, true
647
648 );
649 };
650
651 auto set_ops = [&](auto schur_dm) {
653 auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
654 auto ao_up = createAOMappingIS(dm_is, PETSC_NULLPTR);
655
656 boost::shared_ptr<BlockStructure> block_data;
657 CHKERR DMMoFEMGetBlocMatData(simple->getDM(), block_data);
658
659 if (AT == BLOCK_SCHUR) {
660 pip_mng->getOpDomainLhsPipeline().push_front(
662 pip_mng->getOpDomainLhsPipeline().push_back(
663
664 createOpSchurAssembleEnd({"BROKEN", "U"}, {nullptr, nullptr}, ao_up,
665 S, true, true)
666
667 );
668 }
669
670 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
671 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
672
673 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
675 CHKERR MatZeroEntries(S);
676 MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble Begin";
678 };
679
680 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
681 post_proc_schur_lhs_ptr]() {
683 MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble End";
684
685 if (AT == BLOCK_SCHUR) {
686 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
687 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
688 if (1) {
689 auto S_from_block = matDuplicate(S, MAT_SHARE_NONZERO_PATTERN);
690 // Create matrix from block mat
691 CHKERR assembleBlockMatSchur(mField, post_proc_schur_lhs_ptr->B,
692 S_from_block, {"BROKEN", "U"},
693 {nullptr, nullptr}, ao_up);
694 CHKERR MatAssemblyBegin(S_from_block, MAT_FINAL_ASSEMBLY);
695 CHKERR MatAssemblyEnd(S_from_block, MAT_FINAL_ASSEMBLY);
696 CHKERR MatAYPX(S_from_block, -1, S, DIFFERENT_NONZERO_PATTERN);
697 double norm;
698 CHKERR MatNorm(S_from_block, NORM_FROBENIUS, &norm);
699 MOFEM_LOG("AT", Sev::inform) << "Norm of difference: " << norm;
700 if (norm > 1e-6)
701 SETERRQ(
702 PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
703 "Norm of difference between Schur and block matrix is larger "
704 "than accepted");
705 }
706 } else {
707 CHKERR assembleBlockMatSchur(mField, post_proc_schur_lhs_ptr->B, S,
708 {"BROKEN", "U"}, {nullptr, nullptr}, ao_up);
709 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
710 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
711 }
712
713 MOFEM_LOG("AT", Sev::verbose) << "Lhs Assemble Finish";
715 };
716
717 auto ksp_ctx_ptr = getDMKspCtx(simple->getDM());
718 ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
719 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
720
722 };
723
724 auto set_pc = [&](auto pc, auto block_dm) {
726 auto block_is = getDMSubData(block_dm)->getSmartRowIs();
727 CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
728 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
730 };
731
732
733
734 auto [schur_dm, block_dm] = create_sub_dm();
735 if (AT == BLOCK_SCHUR || AT == BLOCK_MAT) {
736 auto nested_mat_data = get_nested_mat_data(schur_dm, block_dm);
737 CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
738 }
739 S = createDMHybridisedL2Matrix(schur_dm);
740 CHKERR MatSetDM(S, PETSC_NULLPTR);
741
742 int bs = (SPACE_DIM == 2) ? NBEDGE_L2(approx_order - 1)
744 CHKERR MatSetBlockSize(S, bs);
745
746 CHKERR set_ops(schur_dm);
747 CHKERR set_pc(pc, block_dm);
748 DM solver_dm;
749 CHKERR KSPGetDM(ksp, &solver_dm);
750 if (AT == BLOCK_SCHUR || AT == BLOCK_MAT)
751 CHKERR DMSetMatType(solver_dm, MATSHELL);
752
753 auto get_pc = [](auto ksp) {
754 PC pc_raw;
755 CHKERR KSPGetPC(ksp, &pc_raw);
756 return pc_raw;
757 };
758
759 auto set_diagonal_pc = [&](auto pc) {
761
762 if (AT == BLOCK_SCHUR || AT == BLOCK_MAT) {
763 auto A = createDMBlockMat(simple->getDM());
764 auto P = createDMNestSchurMat(simple->getDM());
765 CHKERR PCSetOperators(pc, A, P);
766 }
767
768 KSP *subksp;
769 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
770 CHKERR setSchurA00MatSolvePC(SmartPetscObj<PC>(get_pc(subksp[0]), true));
771
772 CHKERR PetscFree(subksp);
774 };
775
776 auto set_mg_for_schur_complement = [&](auto pc, auto schur_dm, auto S) {
778
779 KSP *subksp;
780 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
781 auto subpc = get_pc(subksp[1]);
782
783 auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
785 CHKERR PCSetDM(pc, dm);
786 PetscBool same = PETSC_FALSE;
787 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
788 if (same) {
789 MOFEM_LOG("TIMER", Sev::inform) << "Set up MG";
791 pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
792 CHKERR PCSetFromOptions(pc);
793 }
795 };
796
797 PetscBool same = PETSC_FALSE;
798 PetscObjectTypeCompare((PetscObject)subpc, PCKSP, &same);
799 if (same) {
800 CHKERR PCSetFromOptions(subpc);
801 KSP inner_ksp;
802 CHKERR PCKSPGetKSP(subpc, &inner_ksp);
803 CHKERR KSPSetFromOptions(inner_ksp);
804 PC inner_pc;
805 CHKERR KSPGetPC(inner_ksp, &inner_pc);
806 CHKERR PCSetFromOptions(inner_pc);
807 CHKERR set_pc_p_mg(schur_dm, inner_pc, S);
808 }
809
810 CHKERR PetscFree(subksp);
812 };
813
814 CHKERR KSPSetUp(ksp);
815 if (AT == BLOCK_SCHUR || AT == BLOCK_MAT) {
816 CHKERR set_diagonal_pc(pc);
817 CHKERR set_mg_for_schur_complement(pc, schur_dm, S);
818 }
819
820 } else {
821 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
822 "PC is not set to PCFIELDSPLIT");
823 }
825}
826
827boost::shared_ptr<SetUpSchur>
829 return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
830}
Integrator for broken space constraints.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
FieldSpace
approximation spaces
Definition definitions.h:82
@ L2
field with C-1 continuity
Definition definitions.h:88
@ NOSPACE
Definition definitions.h:83
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto integration_rule
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
#define NBEDGE_L2(P)
Number of base functions on edge from L2 space.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition Schur.cpp:2585
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode DMMoFEMGetBlocMatData(DM dm, boost::shared_ptr< BlockStructure > &)
Get data for block mat.
Definition DMMoFEM.cpp:1534
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition Schur.cpp:2627
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition DMMoFEM.hpp:1116
auto createDMHybridisedL2Matrix(DM dm)
Get smart hybridised L2 matrix from DM.
Definition DMMoFEM.hpp:1072
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition DMMoFEM.hpp:1160
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition Schur.cpp:1082
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition Schur.cpp:2343
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition DMMoFEM.cpp:1554
auto createDMNestSchurMat(DM dm)
Definition DMMoFEM.hpp:1086
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao)
Assemble Schur matrix.
Definition Schur.cpp:1817
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition Schur.cpp:2580
auto createDMBlockMat(DM dm)
Definition DMMoFEM.hpp:1079
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition seepage.cpp:107
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
Definition seepage.cpp:85
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition seepage.cpp:114
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
MoFEMErrorCode removeSideDOFs(const std::string problem_name, const std::string block_name, const std::string field_name, int bridge_dim, int lo, int hi, bool is_distributed_mesh=true)
Remove side DOFs.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Operator for broken loop side.
Get norm of input VectorDouble for Tensor0.
Get norm of input MatrixDouble for Tensor1.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
virtual ~SetUpSchurImpl()=default
SmartPetscObj< Mat > S
SetUpSchurImpl(MoFEM::Interface &m_field)
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
MoFEM::Interface & mField
[Push operators to pipeline]
SetUpSchur()=default
virtual MoFEMErrorCode setUp(SmartPetscObj< KSP >)=0
virtual ~SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
constexpr AssemblyType AT
int approx_order
constexpr IntegrationType IT
static char help[]
constexpr AssemblyType BlockAssemblyType
constexpr int SPACE_DIM
constexpr bool debug
BoundaryEle::UserDataOperator BdyEleOp