v0.14.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
test_broken_space.cpp File Reference
#include <MoFEM.hpp>
#include <FormsBrokenSpaceConstraintImpl.hpp>

Go to the source code of this file.

Classes

struct  SetUpSchur
 [Push operators to pipeline] More...
 
struct  SetUpSchurImpl
 

Typedefs

using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle
 
using BoundaryEle
 
using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle
 
using BdyEleOp = BoundaryEle::UserDataOperator
 
using SideEleOp = EleOnSide::UserDataOperator
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
constexpr bool debug = false
 
constexpr AssemblyType AT
 
constexpr IntegrationType IT
 
constexpr int SPACE_DIM
 
int approx_order = 1
 

Typedef Documentation

◆ BdyEleOp

using BdyEleOp = BoundaryEle::UserDataOperator
Examples
test_broken_space.cpp.

Definition at line 37 of file test_broken_space.cpp.

◆ BoundaryEle

using BoundaryEle
Initial value:
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle

Definition at line 31 of file test_broken_space.cpp.

◆ DomainEle

using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle

Definition at line 30 of file test_broken_space.cpp.

◆ EleOnSide

using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle

Definition at line 33 of file test_broken_space.cpp.

◆ SideEleOp

using SideEleOp = EleOnSide::UserDataOperator

Definition at line 38 of file test_broken_space.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

[Register MoFEM discrete manager in PETSc]

[Register MoFEM discrete manager in PETSc

Definition at line 52 of file test_broken_space.cpp.

52 {
53
54 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
55
56 try {
57
58 //! [Register MoFEM discrete manager in PETSc]
59 DMType dm_name = "DMMOFEM";
60 CHKERR DMRegister_MoFEM(dm_name);
61 DMType dm_name_mg = "DMMOFEM_MG";
62 CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
63 //! [Register MoFEM discrete manager in PETSc
64
65 moab::Core mb_instance;
66 moab::Interface &moab = mb_instance;
67
68 // Add logging channel for example
69 auto core_log = logging::core::get();
70 core_log->add_sink(
71 LogManager::createSink(LogManager::getStrmWorld(), "AT"));
72 core_log->add_sink(
73 LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
74 LogManager::setLog("AT");
75 LogManager::setLog("TIMER");
76 MOFEM_LOG_TAG("AT", "atom_test");
77 MOFEM_LOG_TAG("TIMER", "timer");
78
79 // Create MoFEM instance
80 MoFEM::Core core(moab);
81 MoFEM::Interface &m_field = core;
82
83 auto *simple = m_field.getInterface<Simple>();
84 CHKERR simple->getOptions();
85 simple->getAddBoundaryFE() = true;
86 CHKERR simple->loadFile();
87
88 auto add_shared_entities_on_skeleton = [&]() {
89 MoFEMFunctionBegin;
90 auto boundary_meshset = simple->getBoundaryMeshSet();
91 auto skeleton_meshset = simple->getSkeletonMeshSet();
92 Range bdy_ents;
93 CHKERR m_field.get_moab().get_entities_by_handle(boundary_meshset,
94 bdy_ents, true);
95 Range skeleton_ents;
96 CHKERR m_field.get_moab().get_entities_by_dimension(
97 0, simple->getDim() - 1, skeleton_ents, true);
98 skeleton_ents = subtract(skeleton_ents, bdy_ents);
99 CHKERR m_field.get_moab().clear_meshset(&skeleton_meshset, 1);
100 CHKERR m_field.get_moab().add_entities(skeleton_meshset, skeleton_ents);
101 MoFEMFunctionReturn(0);
102 };
103
104 CHKERR add_shared_entities_on_skeleton();
105
106 // Declare elements
107 enum bases {
108 AINSWORTH,
109 AINSWORTH_LOBATTO,
110 DEMKOWICZ,
111 BERNSTEIN,
112 LASBASETOP
113 };
114 const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
115 "bernstein"};
116 PetscBool flg;
117 PetscInt choice_base_value = AINSWORTH;
118 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
119 LASBASETOP, &choice_base_value, &flg);
120
121 if (flg != PETSC_TRUE)
122 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
123 FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
124 if (choice_base_value == AINSWORTH)
125 base = AINSWORTH_LEGENDRE_BASE;
126 if (choice_base_value == AINSWORTH_LOBATTO)
127 base = AINSWORTH_LOBATTO_BASE;
128 else if (choice_base_value == DEMKOWICZ)
129 base = DEMKOWICZ_JACOBI_BASE;
130 else if (choice_base_value == BERNSTEIN)
131 base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
132
133 enum spaces { hdiv, hcurl, last_space };
134 const char *list_spaces[] = {"hdiv", "hcurl"};
135 PetscInt choice_space_value = hdiv;
136 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
137 last_space, &choice_space_value, &flg);
138 if (flg != PETSC_TRUE)
139 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
140 FieldSpace space = HDIV;
141 if (choice_space_value == hdiv)
142 space = HDIV;
143 else if (choice_space_value == hcurl)
144 space = HCURL;
145
146 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
147 PETSC_NULL);
148
149 CHKERR simple->addDomainBrokenField("BROKEN", space, base, 1);
150 CHKERR simple->addDomainField("U", L2, base, 1);
151 CHKERR simple->addSkeletonField("HYBRID", L2, base, 1);
152
153 CHKERR simple->setFieldOrder("BROKEN", approx_order);
154 CHKERR simple->setFieldOrder("U", approx_order - 1);
155 CHKERR simple->setFieldOrder("HYBRID", approx_order - 1);
156
157 CHKERR simple->setUp();
158
159 auto bc_mng = m_field.getInterface<BcManager>();
160 CHKERR bc_mng->removeSideDOFs(simple->getProblemName(), "ZERO_FLUX",
161 "BROKEN", SPACE_DIM, 0, 1, true);
162
163 auto integration_rule = [](int, int, int p) { return 2 * p; };
164
165 auto assemble_domain_lhs = [&](auto &pip) {
166 MoFEMFunctionBegin;
167 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
168
169 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<
170 AT>::BiLinearForm<IT>::OpMass<3, SPACE_DIM>;
171 using OpHdivU = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
172 IT>::OpMixDivTimesScalar<SPACE_DIM>;
173
174 auto beta = [](const double, const double, const double) constexpr {
175 return 1;
176 };
177
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));
181
182 // First: Iterate over skeleton FEs adjacent to Domain FEs
183 // Note: BoundaryEle, i.e. uses skeleton interation rule
184 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
185 m_field, simple->getSkeletonFEName(), SPACE_DIM - 1, Sev::noisy);
186 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = integration_rule;
187 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
188 op_loop_skeleton_side->getOpPtrVector(), {});
189
190 // Second: Iterate over domain FEs adjacent to skelton, particularly one
191 // domain element.
192 auto broken_data_ptr =
193 boost::make_shared<std::vector<BrokenBaseSideData>>();
194 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
195 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
196 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy);
197 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
198 op_loop_domain_side->getOpPtrVector(), {HDIV});
199 op_loop_domain_side->getOpPtrVector().push_back(
200 new OpGetBrokenBaseSideData<SideEleOp>("BROKEN", broken_data_ptr));
201
202 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
203 using OpC = FormsIntegrators<BdyEleOp>::Assembly<AT>::BiLinearForm<
204 IT>::OpBrokenSpaceConstrain<1>;
205 op_loop_skeleton_side->getOpPtrVector().push_back(
206 new OpC("HYBRID", broken_data_ptr, 1., true, false));
207
208 if (debug) {
209 // print skeleton elements on partition
210 constexpr int partition = 1;
211 auto op_print = new BdyEleOp(NOSPACE, BdyEleOp::OPSPACE);
212 op_print->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
213 EntityType type,
214 EntitiesFieldData::EntData &data) {
215 MoFEMFunctionBegin;
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;
219
220 if (m_field.get_comm_rank() == partition) {
221 if (num_fe->getPStatus() & PSTATUS_SHARED)
222 MOFEM_LOG("SELF", Sev::inform) << "Num FE: " << *num_fe;
223 }
224 }
225 MoFEMFunctionReturn(0);
226 };
227 op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
228 };
229
230 pip.push_back(op_loop_skeleton_side);
231
232 MoFEMFunctionReturn(0);
233 };
234
235 auto assemble_domain_rhs = [&](auto &pip) {
236 MoFEMFunctionBegin;
237 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
238 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
239 AT>::LinearForm<IT>::OpSource<1, 1>;
240 auto source = [&](const double x, const double y,
241 const double z) constexpr {
242 return -1; // sin(100 * (x / 10.) * M_PI_2);
243 };
244 pip.push_back(new OpDomainSource("U", source));
245 MoFEMFunctionReturn(0);
246 };
247
248 auto *pip_mng = m_field.getInterface<PipelineManager>();
249
250 CHKERR assemble_domain_lhs(pip_mng->getOpDomainLhsPipeline());
251 CHKERR assemble_domain_rhs(pip_mng->getOpDomainRhsPipeline());
252
253 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
254 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
255 CHKERR pip_mng->setSkeletonLhsIntegrationRule(integration_rule);
256 CHKERR pip_mng->setSkeletonRhsIntegrationRule(integration_rule);
257
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()});
262
263 auto x = createDMVector(simple->getDM());
264 auto f = vectorDuplicate(x);
265
266 if (AT == PETSC) {
267 auto ksp = pip_mng->createKSP();
268
269 CHKERR KSPSetFromOptions(ksp);
270 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
271 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
272 CHKERR KSPSetUp(ksp);
273 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
274
275 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
276 CHKERR KSPSolve(ksp, f, x);
277 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
278
279 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
280 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
281 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
282 SCATTER_REVERSE);
283 } else {
284 auto ksp = pip_mng->createKSP();
285 auto schur_ptr = SetUpSchur::createSetUpSchur(m_field);
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";
290
291 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
292 CHKERR KSPSolve(ksp, f, x);
293 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
294
295 CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
296 CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
297 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), x, INSERT_VALUES,
298 SCATTER_REVERSE);
299 }
300
301 auto check_residual = [&](auto x, auto f) {
302 MoFEMFunctionBegin;
303 auto *simple = m_field.getInterface<Simple>();
304 auto *pip_mng = m_field.getInterface<PipelineManager>();
305
306 // auto &skeleton_rhs = pip_mng->getOpSkeletonRhsPipeline();
307 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
308 // skeleton_rhs.clear();
309 domain_rhs.clear();
310
311 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(domain_rhs, {HDIV});
312
313 auto div_flux_ptr = boost::make_shared<VectorDouble>();
314 domain_rhs.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
315 "BROKEN", div_flux_ptr));
316 using OpUDivFlux = FormsIntegrators<DomainEleOp>::Assembly<
317 AT>::LinearForm<IT>::OpBaseTimesScalarField<1>;
318 auto beta = [](double, double, double) constexpr { return 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; };
322 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
323 AT>::LinearForm<IT>::OpSource<1, 1>;
324 domain_rhs.push_back(new OpDomainSource("U", source));
325
326 using OpHDivH = FormsIntegrators<DomainEleOp>::Assembly<AT>::LinearForm<
327 IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
328 using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<
329 AT>::LinearForm<IT>::OpBaseTimesVector<3, 3, 1>;
330 auto flux_ptr = boost::make_shared<MatrixDouble>();
331 domain_rhs.push_back(
332 new OpCalculateHVecVectorField<3>("BROKEN", flux_ptr));
333 boost::shared_ptr<VectorDouble> u_ptr =
334 boost::make_shared<VectorDouble>();
335 domain_rhs.push_back(new OpCalculateScalarFieldValues("U", u_ptr));
336 domain_rhs.push_back(new OpHDivH("BROKEN", u_ptr, beta));
337 domain_rhs.push_back(new OpHdivFlux("BROKEN", flux_ptr, beta));
338
339 // First: Iterate over skeleton FEs adjacent to Domain FEs
340 // Note: BoundaryEle, i.e. uses skeleton interation rule
341 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
342 m_field, simple->getSkeletonFEName(), SPACE_DIM - 1, Sev::noisy);
343 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = integration_rule;
344 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
345 op_loop_skeleton_side->getOpPtrVector(), {});
346
347 // Second: Iterate over domain FEs adjacent to skelton, particularly one
348 // domain element.
349 auto broken_data_ptr =
350 boost::make_shared<std::vector<BrokenBaseSideData>>();
351 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
352 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
353 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy);
354 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
355 op_loop_domain_side->getOpPtrVector(), {HDIV});
356 op_loop_domain_side->getOpPtrVector().push_back(
357 new OpGetBrokenBaseSideData<SideEleOp>("BROKEN", broken_data_ptr));
358 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
359 op_loop_domain_side->getOpPtrVector().push_back(
360 new OpCalculateHVecTensorField<1, 3>("BROKEN", flux_mat_ptr));
361 op_loop_domain_side->getOpPtrVector().push_back(
362 new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
363
364 // Assemble on skeleton
365 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
366 using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<AT>::LinearForm<
367 IT>::OpBrokenSpaceConstrainDHybrid<1>;
368 using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<AT>::LinearForm<
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(
374 new OpCalculateVectorFieldValues<1>("HYBRID", hybrid_ptr));
375 op_loop_skeleton_side->getOpPtrVector().push_back(
376 new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
377
378 // Add skeleton to domain pipeline
379 domain_rhs.push_back(op_loop_skeleton_side);
380
381 CHKERR VecZeroEntries(f);
382 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
383 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
384
385 pip_mng->getDomainRhsFE()->f = f;
386 pip_mng->getSkeletonRhsFE()->f = f;
387 pip_mng->getDomainRhsFE()->x = x;
388 pip_mng->getSkeletonRhsFE()->x = x;
389
390 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
391 simple->getDomainFEName(),
392 pip_mng->getDomainRhsFE());
393
394 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
395 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
396 CHKERR VecAssemblyBegin(f);
397 CHKERR VecAssemblyEnd(f);
398
399 double fnrm;
400 CHKERR VecNorm(f, NORM_2, &fnrm);
401 MOFEM_LOG_C("AT", Sev::inform, "Residual %3.4e", fnrm);
402
403 constexpr double eps = 1e-8;
404 if (fnrm > eps)
405 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
406 "Residual norm larger than accepted");
407
408 MoFEMFunctionReturn(0);
409 };
410
411 auto calculate_error = [&]() {
412 MoFEMFunctionBegin;
413
414 // auto &skeleton_rhs = pip_mng->getOpSkeletonRhsPipeline();
415 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
416 // skeleton_rhs.clear();
417 domain_rhs.clear();
418
419 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(domain_rhs, {HDIV});
420
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>();
425
426 domain_rhs.push_back(
427 new OpCalculateScalarFieldGradient<SPACE_DIM>("U", u_grad_ptr));
428 domain_rhs.push_back(
429 new OpCalculateHVecVectorField<3, SPACE_DIM>("BROKEN", flux_val_ptr));
430 domain_rhs.push_back(new OpCalculateHdivVectorDivergence<3, SPACE_DIM>(
431 "BROKEN", div_val_ptr));
432 auto source = [&](const double x, const double y,
433 const double z) constexpr { return -1; };
434 domain_rhs.push_back(new OpGetTensor0fromFunc(source_ptr, source));
435
436 enum { DIV, GRAD, LAST };
437 auto mpi_vec = createVectorMPI(
438 m_field.get_comm(), (!m_field.get_comm_rank()) ? LAST : 0, LAST);
439 domain_rhs.push_back(
440 new OpCalcNormL2Tensor0(div_val_ptr, mpi_vec, DIV, source_ptr));
441 domain_rhs.push_back(new OpCalcNormL2Tensor1<SPACE_DIM>(
442 u_grad_ptr, mpi_vec, GRAD, flux_val_ptr));
443
444 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
445 simple->getDomainFEName(),
446 pip_mng->getDomainRhsFE());
447 CHKERR VecAssemblyBegin(mpi_vec);
448 CHKERR VecAssemblyEnd(mpi_vec);
449
450 if (!m_field.get_comm_rank()) {
451 const double *error_ind;
452 CHKERR VecGetArrayRead(mpi_vec, &error_ind);
453 MOFEM_LOG("AT", Sev::inform)
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);
459 }
460
461 MoFEMFunctionReturn(0);
462 };
463
464 auto get_post_proc_fe = [&]() {
465 using PostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
466 using OpPPMap = OpPostProcMapInMoab<3, SPACE_DIM>;
467 auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
468
469 auto op_loop_side = new OpLoopSide<EleOnSide>(
470 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy,
471 boost::make_shared<
472 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
473 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
474
475 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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(
480 new OpCalculateScalarFieldValues("U", u_vec_ptr));
481 op_loop_side->getOpPtrVector().push_back(
482 new OpCalculateHVecVectorField<3>("BROKEN", flux_mat_ptr));
483 op_loop_side->getOpPtrVector().push_back(
484
485 new OpPPMap(
486
487 post_proc_fe->getPostProcMesh(),
488
489 post_proc_fe->getMapGaussPts(),
490
491 {{"U", u_vec_ptr}},
492
493 {{"BROKEN", flux_mat_ptr}},
494
495 {}, {})
496
497 );
498
499 return post_proc_fe;
500 };
501
502 auto post_proc_fe = get_post_proc_fe();
503 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
504 simple->getBoundaryFEName(), post_proc_fe);
505 CHKERR post_proc_fe->writeFile("out_result.h5m");
506
507 CHKERR calculate_error();
508 CHKERR check_residual(x, f);
509 }
510 CATCH_ERRORS;
511
512 CHKERR MoFEM::Core::Finalize();
513}
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition child_and_parent.cpp:54
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
FieldSpace
FieldSpace
approximation spaces
Definition definitions.h:82
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
NOSPACE
@ NOSPACE
Definition definitions.h:83
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:429
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:548
integration_rule
auto integration_rule
Definition free_surface.cpp:185
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition PCMGSetUpViaApproxOrders.cpp:302
MoFEM::DMoFEMLoopFiniteElements
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:586
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1099
MoFEM::PETSC
@ PETSC
Definition FormsIntegrators.hpp:105
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:389
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:308
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:339
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition mixed_poisson.cpp:24
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:221
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition PetscSmartObj.hpp:202
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition DeprecatedPetsc.hpp:203
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:29
source
auto source
Definition poisson_2d_dis_galerkin.cpp:62
OpHDivH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition seepage.cpp:115
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition HODataOperators.hpp:415
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
MoFEM::BcManager::removeSideDOFs
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.
Definition BcManager.cpp:1835
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition DataOperators.hpp:24
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:128
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition LogManager.cpp:298
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:344
MoFEM::OpCalcNormL2Tensor0
Get norm of input VectorDouble for Tensor0.
Definition NormsOperators.hpp:15
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition NormsOperators.hpp:44
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition UserDataOperators.hpp:2573
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition UserDataOperators.hpp:2117
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition UserDataOperators.hpp:2213
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1295
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition UserDataOperators.hpp:83
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition UserDataOperators.hpp:468
MoFEM::OpGetTensor0fromFunc
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
Definition NormsOperators.hpp:105
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition ForcesAndSourcesCore.hpp:1291
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PipelineManager
PipelineManager interface.
Definition PipelineManager.hpp:24
MoFEM::PostProcBrokenMeshInMoab
Definition PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition test_broken_space.cpp:768
AT
constexpr AssemblyType AT
Definition test_broken_space.cpp:20
approx_order
int approx_order
Definition test_broken_space.cpp:50
IT
constexpr IntegrationType IT
Definition test_broken_space.cpp:24
help
static char help[]
Definition test_broken_space.cpp:16
SPACE_DIM
constexpr int SPACE_DIM
Definition test_broken_space.cpp:27
debug
constexpr bool debug
Definition test_broken_space.cpp:18
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition test_broken_space.cpp:37
OpHdivFlux
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, SPACE_DIM, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition thermo_elastic.cpp:72
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition thermo_elastic.cpp:50

Variable Documentation

◆ approx_order

int approx_order = 1

Definition at line 50 of file test_broken_space.cpp.

◆ AT

AssemblyType AT
constexpr
Initial value:
=
(SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
: AssemblyType::PETSC
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition FormsIntegrators.hpp:104
Examples
plastic.cpp, and test_broken_space.cpp.

Definition at line 20 of file test_broken_space.cpp.

◆ debug

bool debug = false
constexpr

Definition at line 18 of file test_broken_space.cpp.

◆ help

char help[] = "...\n\n"
static

Definition at line 16 of file test_broken_space.cpp.

◆ IT

IntegrationType IT
constexpr
Initial value:
=
IntegrationType::GAUSS
Examples
plastic.cpp, and test_broken_space.cpp.

Definition at line 24 of file test_broken_space.cpp.

◆ SPACE_DIM

int SPACE_DIM
constexpr
Initial value:
=
EXECUTABLE_DIMENSION
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13

Definition at line 27 of file test_broken_space.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow