56 {
57
59
60 try {
61
62
63 DMType dm_name = "DMMOFEM";
65 DMType dm_name_mg = "DMMOFEM_MG";
67
68
69 moab::Core mb_instance;
70 moab::Interface &moab = mb_instance;
71
72
73 auto core_log = logging::core::get();
74 core_log->add_sink(
76 core_log->add_sink(
82
83
86
89 simple->getAddBoundaryFE() =
true;
91
92 auto add_shared_entities_on_skeleton = [&]() {
94 auto boundary_meshset =
simple->getBoundaryMeshSet();
95 auto skeleton_meshset =
simple->getSkeletonMeshSet();
98 bdy_ents, true);
101 0,
simple->getDim() - 1, skeleton_ents,
true);
102 skeleton_ents = subtract(skeleton_ents, bdy_ents);
104 CHKERR m_field.
get_moab().add_entities(skeleton_meshset, skeleton_ents);
106 };
107
108 CHKERR add_shared_entities_on_skeleton();
109
110
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;
123 LASBASETOP, &choice_base_value, &flg);
124
125 if (flg != PETSC_TRUE)
128 if (choice_base_value == AINSWORTH)
130 if (choice_base_value == AINSWORTH_LOBATTO)
132 else if (choice_base_value == DEMKOWICZ)
134 else if (choice_base_value == BERNSTEIN)
136
137 enum spaces { hdiv, hcurl, last_space };
138 const char *list_spaces[] = {"hdiv", "hcurl"};
139 PetscInt choice_space_value = hdiv;
141 last_space, &choice_space_value, &flg);
142 if (flg != PETSC_TRUE)
145 if (choice_space_value == hdiv)
147 else if (choice_space_value == hcurl)
149
151 PETSC_NULLPTR);
152
153 CHKERR simple->addDomainBrokenField(
"BROKEN", space, base, 1);
156
160
162
166
168
169 auto assemble_domain_lhs = [&](auto &pip) {
172
176 IT>::OpMixDivTimesScalar<SPACE_DIM>;
177
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
187
192 op_loop_skeleton_side->getOpPtrVector(), {});
193
194
195
196 auto broken_data_ptr =
197 boost::make_shared<std::vector<BrokenBaseSideData>>();
198
202 op_loop_domain_side->getOpPtrVector(), {HDIV});
203 op_loop_domain_side->getOpPtrVector().push_back(
205
206 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
208 IT>::OpBrokenSpaceConstrain<1>;
209 op_loop_skeleton_side->getOpPtrVector().push_back(
210 new OpC("HYBRID", broken_data_ptr, 1., true, false));
211
213
214 constexpr int partition = 1;
216 op_print->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
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
225 if (num_fe->getPStatus() & PSTATUS_SHARED)
226 MOFEM_LOG(
"SELF", Sev::inform) <<
"Num FE: " << *num_fe;
227 }
228 }
230 };
231 op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
232 };
233
234 pip.push_back(op_loop_skeleton_side);
235
237 };
238
239 auto assemble_domain_rhs = [&](auto &pip) {
244 auto source = [&](
const double x,
const double y,
245 const double z) constexpr {
246 return -1;
247 };
250 };
251
253
254 CHKERR assemble_domain_lhs(pip_mng->getOpDomainLhsPipeline());
255 CHKERR assemble_domain_rhs(pip_mng->getOpDomainRhsPipeline());
256
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
269
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";
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);
286 SCATTER_REVERSE);
287 } else {
288 auto ksp = pip_mng->createKSP();
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);
302 SCATTER_REVERSE);
303 }
304
305 auto check_residual = [&](
auto x,
auto f) {
309
310
311 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
312
313 domain_rhs.clear();
314
316
317 auto div_flux_ptr = boost::make_shared<VectorDouble>();
319 "BROKEN", div_flux_ptr));
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; };
329
331 IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
334 auto flux_ptr = boost::make_shared<MatrixDouble>();
335 domain_rhs.push_back(
337 boost::shared_ptr<VectorDouble> u_ptr =
338 boost::make_shared<VectorDouble>();
340 domain_rhs.push_back(
new OpHDivH(
"BROKEN", u_ptr, beta));
341 domain_rhs.push_back(
new OpHdivFlux(
"BROKEN", flux_ptr, beta));
342
343
344
349 op_loop_skeleton_side->getOpPtrVector(), {});
350
351
352
353 auto broken_data_ptr =
354 boost::make_shared<std::vector<BrokenBaseSideData>>();
355
359 op_loop_domain_side->getOpPtrVector(), {HDIV});
360 op_loop_domain_side->getOpPtrVector().push_back(
362 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
363 op_loop_domain_side->getOpPtrVector().push_back(
365 op_loop_domain_side->getOpPtrVector().push_back(
367
368
369 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
371 IT>::OpBrokenSpaceConstrainDHybrid<1>;
373 IT>::OpBrokenSpaceConstrainDFlux<1>;
374 op_loop_skeleton_side->getOpPtrVector().push_back(
375 new OpC_dHybrid("HYBRID", broken_data_ptr, 1.));
376 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
377 op_loop_skeleton_side->getOpPtrVector().push_back(
379 op_loop_skeleton_side->getOpPtrVector().push_back(
380 new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
381
382
383 domain_rhs.push_back(op_loop_skeleton_side);
384
386 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
387 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
388
389 pip_mng->getDomainRhsFE()->f =
f;
390 pip_mng->getSkeletonRhsFE()->f =
f;
391 pip_mng->getDomainRhsFE()->x = x;
392 pip_mng->getSkeletonRhsFE()->x = x;
393
395 simple->getDomainFEName(),
396 pip_mng->getDomainRhsFE());
397
398 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
399 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
400 CHKERR VecAssemblyBegin(f);
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;
410 "Residual norm larger than accepted");
411
413 };
414
415 auto calculate_error = [&]() {
417
418
419 auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
420
421 domain_rhs.clear();
422
424
425 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
426 auto flux_val_ptr = boost::make_shared<MatrixDouble>();
427 auto div_val_ptr = boost::make_shared<VectorDouble>();
428 auto source_ptr = boost::make_shared<VectorDouble>();
429
430 domain_rhs.push_back(
432 domain_rhs.push_back(
435 "BROKEN", div_val_ptr));
436 auto source = [&](
const double x,
const double y,
437 const double z) constexpr { return -1; };
439
440 enum { DIV, GRAD, LAST };
443 domain_rhs.push_back(
446 u_grad_ptr, mpi_vec, GRAD, flux_val_ptr));
447
449 simple->getDomainFEName(),
450 pip_mng->getDomainRhsFE());
451 CHKERR VecAssemblyBegin(mpi_vec);
452 CHKERR VecAssemblyEnd(mpi_vec);
453
455 const double *error_ind;
456 CHKERR VecGetArrayRead(mpi_vec, &error_ind);
458 << "Approximation error ||div flux - source||: "
459 << std::sqrt(error_ind[DIV]);
460 MOFEM_LOG(
"AT", Sev::inform) <<
"Approximation error ||grad-flux||: "
461 << std::sqrt(error_ind[GRAD]);
462 CHKERR VecRestoreArrayRead(mpi_vec, &error_ind);
463 }
464
466 };
467
468 auto get_post_proc_fe = [&]() {
471 auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
472
475 boost::make_shared<
476 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
477 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
478
480 op_loop_side->getOpPtrVector(), {HDIV});
481 auto u_vec_ptr = boost::make_shared<VectorDouble>();
482 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
483 op_loop_side->getOpPtrVector().push_back(
485 op_loop_side->getOpPtrVector().push_back(
487 op_loop_side->getOpPtrVector().push_back(
488
490
491 post_proc_fe->getPostProcMesh(),
492
493 post_proc_fe->getMapGaussPts(),
494
495 {{"U", u_vec_ptr}},
496
497 {{"BROKEN", flux_mat_ptr}},
498
499 {}, {})
500
501 );
502
503 return post_proc_fe;
504 };
505
506 auto post_proc_fe = get_post_proc_fe();
508 simple->getBoundaryFEName(), post_proc_fe);
509 CHKERR post_proc_fe->writeFile(
"out_result.h5m");
510
512 CHKERR check_residual(x, f);
513 }
515
517}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, 2 > OpHDivH
Integrate Rhs div flux base times temperature (T)
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
MoFEMErrorCode removeSideDOFs(const std::string problem_name, const std::string block_name, const std::string field_name, int bridge_dim, int lo, int hi, bool is_distributed_mesh=true)
Remove side DOFs.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Operator for broken loop side.
Get norm of input VectorDouble for Tensor0.
Get norm of input MatrixDouble for Tensor1.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
PipelineManager interface.
Simple interface for fast problem set-up.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
constexpr AssemblyType AT
constexpr IntegrationType IT
BoundaryEle::UserDataOperator BdyEleOp