v0.14.0
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
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 36 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
Examples
test_broken_space.cpp.

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

Examples
test_broken_space.cpp.

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_NULL, 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_NULL, 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_NULL, "", "-order", &approx_order,
151  PETSC_NULL);
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 }

Variable Documentation

◆ approx_order

int approx_order = 1
Examples
adolc_plasticity.cpp, approx_sphere.cpp, dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_equation.cpp, photon_diffusion.cpp, prism_polynomial_approximation.cpp, quad_polynomial_approximation.cpp, scalar_check_approximation.cpp, seepage.cpp, shallow_wave.cpp, test_broken_space.cpp, thermo_elastic.cpp, and wave_equation.cpp.

Definition at line 54 of file test_broken_space.cpp.

◆ AT

constexpr AssemblyType AT
constexpr
Initial value:
=
(SCHUR_ASSEMBLE) ? BlockAssemblyType
: AssemblyType::PETSC
Examples
test_broken_space.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
Examples
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
Examples
test_broken_space.cpp.

Definition at line 16 of file test_broken_space.cpp.

◆ IT

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

Definition at line 28 of file test_broken_space.cpp.

◆ SPACE_DIM

constexpr int SPACE_DIM
constexpr
Initial value:
=
EXECUTABLE_DIMENSION
Examples
test_broken_space.cpp.

Definition at line 31 of file test_broken_space.cpp.

NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SPACE_DIM
constexpr int SPACE_DIM
Definition: test_broken_space.cpp:31
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::ForcesAndSourcesCore::UserDataOperator::AdjCache
std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > >> AdjCache
Definition: ForcesAndSourcesCore.hpp:904
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
debug
constexpr bool debug
Definition: test_broken_space.cpp:18
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition: thermo_elastic.cpp:51
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
BlockAssemblyType
constexpr AssemblyType BlockAssemblyType
Definition: test_broken_space.cpp:20
MoFEM::BLOCK_MAT
@ BLOCK_MAT
Definition: FormsIntegrators.hpp:107
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
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::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::OpSetFlux
Definition: FormsBrokenSpaceConstraintImpl.hpp:139
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1319
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
MoFEM::OpGetBrokenBaseSideData
Definition: FormsBrokenSpaceConstraintImpl.hpp:68
OpHDivH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Definition: seepage.cpp:116
AT
constexpr AssemblyType AT
Definition: test_broken_space.cpp:24
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::OpBrokenLoopSide
Definition: FormsBrokenSpaceConstraintImpl.hpp:15
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2171
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
double
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2628
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::OpGetTensor0fromFunc
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
Definition: NormsOperators.hpp:105
BiLinearForm
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
Range
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition: mixed_poisson.cpp:25
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
approx_order
int approx_order
Definition: test_broken_space.cpp:54
IT
constexpr IntegrationType IT
Definition: test_broken_space.cpp:28
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpCalcNormL2Tensor0
Get norm of input VectorDouble for Tensor0.
Definition: NormsOperators.hpp:15
MoFEM::MPC::LAST
@ LAST
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
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
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition: UserDataOperators.hpp:2268
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition: test_broken_space.cpp:41
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
help
static char help[]
Definition: test_broken_space.cpp:16
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
source
auto source
Definition: poisson_2d_dis_galerkin.cpp:62
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:802
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
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
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
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:73
Generated by Doxygen 1.8.17 and hosted at University of Glasgow