v0.15.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 = PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::BoundaryEle
 
using EleOnSide = PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle
 
using EntData = EntitiesFieldData::EntData
 
using DomainEleOp = DomainEle::UserDataOperator
 
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 BlockAssemblyType
 
constexpr AssemblyType AT
 
constexpr IntegrationType IT
 
constexpr int SPACE_DIM
 
int approx_order = 1
 

Typedef Documentation

◆ BdyEleOp

using BdyEleOp = BoundaryEle::UserDataOperator
Examples
mofem/atom_tests/test_broken_space.cpp.

Definition at line 41 of file test_broken_space.cpp.

◆ BoundaryEle

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

Definition at line 35 of file test_broken_space.cpp.

◆ DomainEle

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

Definition at line 34 of file test_broken_space.cpp.

◆ DomainEleOp

using DomainEleOp = DomainEle::UserDataOperator

Definition at line 40 of file test_broken_space.cpp.

◆ EleOnSide

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

Definition at line 37 of file test_broken_space.cpp.

◆ EntData

using EntData = EntitiesFieldData::EntData

Definition at line 39 of file test_broken_space.cpp.

◆ SideEleOp

using SideEleOp = EleOnSide::UserDataOperator

Definition at line 42 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 56 of file test_broken_space.cpp.

56 {
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";
66 CHKERR DMRegister_MGViaApproxOrders(dm_name_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(
75 LogManager::createSink(LogManager::getStrmWorld(), "AT"));
76 core_log->add_sink(
77 LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
78 LogManager::setLog("AT");
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>();
88 CHKERR simple->getOptions();
89 simple->getAddBoundaryFE() = true;
90 CHKERR simple->loadFile();
91
92 auto add_shared_entities_on_skeleton = [&]() {
93 MoFEMFunctionBegin;
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);
105 MoFEMFunctionReturn(0);
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");
127 FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
128 if (choice_base_value == AINSWORTH)
129 base = AINSWORTH_LEGENDRE_BASE;
130 if (choice_base_value == AINSWORTH_LOBATTO)
131 base = AINSWORTH_LOBATTO_BASE;
132 else if (choice_base_value == DEMKOWICZ)
133 base = DEMKOWICZ_JACOBI_BASE;
134 else if (choice_base_value == BERNSTEIN)
135 base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
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) {
170 MoFEMFunctionBegin;
171 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
172
173 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<
174 AT>::BiLinearForm<IT>::OpMass<3, SPACE_DIM>;
175 using OpHdivU = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
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;
191 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
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);
201 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
207 using OpC = FormsIntegrators<BdyEleOp>::Assembly<AT>::BiLinearForm<
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,
218 EntitiesFieldData::EntData &data) {
219 MoFEMFunctionBegin;
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 }
229 MoFEMFunctionReturn(0);
230 };
231 op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
232 };
233
234 pip.push_back(op_loop_skeleton_side);
235
236 MoFEMFunctionReturn(0);
237 };
238
239 auto assemble_domain_rhs = [&](auto &pip) {
240 MoFEMFunctionBegin;
241 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV});
242 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
243 AT>::LinearForm<IT>::OpSource<1, 1>;
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));
249 MoFEMFunctionReturn(0);
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) {
306 MoFEMFunctionBegin;
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
315 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(domain_rhs, {HDIV});
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));
320 using OpUDivFlux = FormsIntegrators<DomainEleOp>::Assembly<
321 AT>::LinearForm<IT>::OpBaseTimesScalarField<1>;
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; };
326 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
327 AT>::LinearForm<IT>::OpSource<1, 1>;
328 domain_rhs.push_back(new OpDomainSource("U", source));
329
330 using OpHDivH = FormsIntegrators<DomainEleOp>::Assembly<AT>::LinearForm<
331 IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
332 using OpHdivFlux = FormsIntegrators<DomainEleOp>::Assembly<
333 AT>::LinearForm<IT>::OpBaseTimesVector<3, 3, 1>;
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;
348 CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
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);
358 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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);
370 using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<AT>::LinearForm<
371 IT>::OpBrokenSpaceConstrainDHybrid<1>;
372 using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<AT>::LinearForm<
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
394 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
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
412 MoFEMFunctionReturn(0);
413 };
414
415 auto calculate_error = [&]() {
416 MoFEMFunctionBegin;
417
418 // auto &skeleton_rhs = pip_mng->getOpSkeletonRhsPipeline();
419 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
420 // skeleton_rhs.clear();
421 domain_rhs.clear();
422
423 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(domain_rhs, {HDIV});
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
448 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
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
465 MoFEMFunctionReturn(0);
466 };
467
468 auto get_post_proc_fe = [&]() {
469 using PostProcEle = PostProcBrokenMeshInMoab<BoundaryEle>;
470 using OpPPMap = OpPostProcMapInMoab<3, SPACE_DIM>;
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
479 CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
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();
507 CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
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 }
514 CATCH_ERRORS;
515
516 CHKERR MoFEM::Core::Finalize();
517}
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
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:56
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 .
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:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
integration_rule
auto integration_rule
Definition free_surface.cpp:211
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:514
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:576
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
MoFEM::PETSC
@ PETSC
Standard PETSc assembly.
Definition FormsIntegrators.hpp:201
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
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 DOFs on side entities from problem.
Definition BcManager.cpp:1686
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition mixed_poisson.cpp:25
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:223
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition PetscSmartObj.hpp:204
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
OpHdivFlux
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
Definition seepage.cpp:108
OpHdivHdiv
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:86
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:599
MoFEM::BcManager
Boundary condition manager for finite element problem setup.
Definition BcManager.hpp:384
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:118
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:152
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:293
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition LogManager.cpp:339
MoFEM::OpBrokenLoopSide
Operator for broken loop side.
Definition FormsBrokenSpaceConstraintImpl.hpp:33
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:3020
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition UserDataOperators.hpp:2402
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition UserDataOperators.hpp:2498
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1490
MoFEM::OpCalculateScalarFieldValues
Specialization for double precision scalar field values calculation.
Definition UserDataOperators.hpp:204
MoFEM::OpCalculateVectorFieldValues
Specialization for double precision vector field values calculation.
Definition UserDataOperators.hpp:626
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:1399
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::PipelineManager
PipelineManager interface.
Definition PipelineManager.hpp:24
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
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:828
AT
constexpr AssemblyType AT
Definition test_broken_space.cpp:24
approx_order
int approx_order
Definition test_broken_space.cpp:54
IT
constexpr IntegrationType IT
Definition test_broken_space.cpp:28
help
static char help[]
Definition test_broken_space.cpp:16
SPACE_DIM
constexpr int SPACE_DIM
Definition test_broken_space.cpp:31
debug
constexpr bool debug
Definition test_broken_space.cpp:18
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition test_broken_space.cpp:41

Variable Documentation

◆ approx_order

int approx_order = 1

Definition at line 54 of file test_broken_space.cpp.

◆ AT

constexpr AssemblyType AT
constexpr
Initial value:
=
(SCHUR_ASSEMBLE) ? BlockAssemblyType
: AssemblyType::PETSC
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition contact.cpp:18
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition FormsIntegrators.hpp:200
BlockAssemblyType
constexpr AssemblyType BlockAssemblyType
Definition test_broken_space.cpp:20
Examples
mofem/atom_tests/test_broken_space.cpp, mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/adv-2/thermo_elastic.cpp.

Definition at line 24 of file test_broken_space.cpp.

◆ BlockAssemblyType

constexpr AssemblyType BlockAssemblyType
constexpr
Initial value:
= (BLOCK_ASSEMBLE_SELECTION)
? AssemblyType::BLOCK_MAT
: AssemblyType::BLOCK_SCHUR
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Block Schur assembly.
Definition FormsIntegrators.hpp:204
Examples
mofem/atom_tests/test_broken_space.cpp.

Definition at line 20 of file test_broken_space.cpp.

◆ debug

constexpr 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

constexpr IntegrationType IT
constexpr
Initial value:
=
IntegrationType::GAUSS
Examples
mofem/atom_tests/test_broken_space.cpp, mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/adv-2/thermo_elastic.cpp.

Definition at line 28 of file test_broken_space.cpp.

◆ SPACE_DIM

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

Definition at line 31 of file test_broken_space.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow