296 {
297
298
301
302 try {
303
304
305 moab::Core mb_instance;
306 moab::Interface &moab = mb_instance;
309
310
311 DMType dm_name = "DMMOFEM";
313
314
319
322
323
325 1);
326
328
330
331
332 boost::shared_ptr<CommonData> data(
new CommonData());
333
334 auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
335 auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
336 auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
337
338
339
340 boost::shared_ptr<Ele> vol_ele_slow_rhs(
new Ele(m_field));
341 boost::shared_ptr<Ele> vol_ele_stiff_rhs(
new Ele(m_field));
342 boost::shared_ptr<Ele> vol_ele_stiff_lhs(
new Ele(m_field));
343
344
345 vol_ele_slow_rhs->getOpPtrVector().push_back(
348
349
350
351
352 auto solve_for_g = [&]() {
354 if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
355 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
356 SCATTER_REVERSE);
357 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
358 SCATTER_REVERSE);
359 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
360 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
361 *vol_ele_slow_rhs->vecAssembleSwitch = false;
362 }
363 CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
364 vol_ele_slow_rhs->ts_F);
366 };
367
368 vol_ele_slow_rhs->postProcessHook = solve_for_g;
369
370 auto det_ptr = boost::make_shared<VectorDouble>();
371 auto jac_ptr = boost::make_shared<MatrixDouble>();
372 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
373
374 vol_ele_stiff_rhs->getOpPtrVector().push_back(
376 vol_ele_stiff_rhs->getOpPtrVector().push_back(
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
380 vol_ele_stiff_rhs->getOpPtrVector().push_back(
382 vol_ele_stiff_rhs->getOpPtrVector().push_back(
384 vol_ele_stiff_rhs->getOpPtrVector().push_back(
386
387
388 vol_ele_stiff_lhs->getOpPtrVector().push_back(
390 vol_ele_stiff_lhs->getOpPtrVector().push_back(
392 vol_ele_stiff_lhs->getOpPtrVector().push_back(
394 vol_ele_stiff_lhs->getOpPtrVector().push_back(
396
397
398 auto vol_rule = [](int, int,
int p) ->
int {
return 2 *
p; };
399 vol_ele_slow_rhs->getRuleHook = vol_rule;
400 vol_ele_stiff_rhs->getRuleHook = vol_rule;
401 vol_ele_stiff_lhs->getRuleHook = vol_rule;
402
403
404
405 auto post_proc = boost::make_shared<PostProcEle>(m_field);
406 boost::shared_ptr<ForcesAndSourcesCore> null;
407
409
410 auto u_ptr = boost::make_shared<VectorDouble>();
411 post_proc->getOpPtrVector().push_back(
413 post_proc->getOpPtrVector().push_back(
414
416
417 post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
418
419 {{"u", u_ptr}},
420
421 {},
422
423 {},
424
425 {}
426
427 )
428
429 );
430
431
432 auto dm = simple_interface->
getDM();
433
434
435
436
440 1,
BLOCKSET, 2, inner_surface,
true);
441 if (!inner_surface.empty()) {
442 Range inner_surface_verts;
443 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts,
false);
445 u0, MBVERTEX, inner_surface_verts,
"u");
446 }
447 }
448
449
450
457 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
458 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
459 PSTATUS_NOT, -1, &edges_part);
461 CHKERR moab.get_connectivity(edges_part, edges_verts,
false);
462
463
466 unite(edges_verts, edges_part));
467
468
470 CHKERR MatZeroEntries(data->M);
471 boost::shared_ptr<Ele> vol_mass_ele(
new Ele(m_field));
472 vol_mass_ele->getOpPtrVector().push_back(
new OpAssembleMass(data));
474 vol_mass_ele);
475 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
476 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
477
478
479
481 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
482 CHKERR KSPSetFromOptions(data->ksp);
483 CHKERR KSPSetUp(data->ksp);
484
485
487
488 CHKERR TSSetType(ts, TSARKIMEX);
489 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
490
491
493 vol_ele_stiff_lhs, null, null);
494
496 vol_ele_stiff_rhs, null, null);
497
499 vol_ele_slow_rhs, null, null);
500
501
502 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, post_proc));
504 monitor_ptr, null, null);
505
506
510
511
512 double ftime = 1;
514 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
515 CHKERR TSSetSolution(ts, X);
516 CHKERR TSSetFromOptions(ts);
518 }
520
521
523
524 return 0;
525}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS the right hand side function
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS(MPI_Comm comm)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
FaceElementForcesAndSourcesCore Ele
const double u0
inital vale on blocksets
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() 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.
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
const std::string getDomainFEName() const
Get the Domain FE Name.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
Assemble stiff part tangent.