2 #include <BasicFiniteElements.hpp>
8 static char help[] =
"...\n\n";
17 , cOmm(m_field.get_comm())
18 , rAnk(m_field.get_comm_rank()) {
19 vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(
new FaceEle(m_field));
20 natural_bdry_ele_slow_rhs =
21 boost::shared_ptr<BoundaryEle>(
new BoundaryEle(m_field));
22 vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(
new FaceEle(m_field));
23 vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(
new FaceEle(m_field));
24 post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
27 data1 = boost::shared_ptr<PreviousData>(
new PreviousData());
28 data2 = boost::shared_ptr<PreviousData>(
new PreviousData());
29 data3 = boost::shared_ptr<PreviousData>(
new PreviousData());
32 boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
34 flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
37 boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
39 mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
42 boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
43 flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
46 boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
47 mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
50 boost::shared_ptr<MatrixDouble>(data3, &data3->flux_values);
51 flux_divs_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->flux_divs);
54 boost::shared_ptr<VectorDouble>(data3, &data3->mass_values);
56 mass_dots_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->mass_dots);
66 MoFEMErrorCode add_fe(std::string mass_field, std::string flux_field);
68 MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL, std::string
internal);
71 boost::shared_ptr<VectorDouble> &mass_ptr);
72 MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
73 boost::shared_ptr<PreviousData> &data);
75 boost::shared_ptr<PreviousData> &data);
78 boost::shared_ptr<VectorDouble> &mass_ptr,
79 boost::shared_ptr<MatrixDouble> &flux_ptr,
80 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
81 boost::shared_ptr<VectorDouble> &flux_div_ptr);
82 MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
83 boost::shared_ptr<PreviousData> &data,
86 std::string flux_field,
87 boost::shared_ptr<VectorDouble> &mass_ptr,
88 boost::shared_ptr<MatrixDouble> &flux_ptr);
89 MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
90 boost::shared_ptr<PreviousData> &data,
95 boost::shared_ptr<FaceEle> &initial_ele);
99 std::string flux_field);
108 Range essential_bdry_ents;
109 Range natural_bdry_ents;
113 Range inner_surface1;
114 Range inner_surface2;
115 Range inner_surface3;
123 std::map<int, BlockData> material_blocks;
125 boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
126 boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
127 boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
128 boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
130 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
131 boost::shared_ptr<Monitor> monitor_ptr;
133 boost::shared_ptr<PreviousData> data1;
134 boost::shared_ptr<PreviousData> data2;
135 boost::shared_ptr<PreviousData>
data3;
137 boost::shared_ptr<MatrixDouble> flux_values_ptr1;
138 boost::shared_ptr<MatrixDouble> flux_values_ptr2;
141 boost::shared_ptr<VectorDouble> flux_divs_ptr1;
142 boost::shared_ptr<VectorDouble> flux_divs_ptr2;
145 boost::shared_ptr<VectorDouble> mass_values_ptr1;
146 boost::shared_ptr<VectorDouble> mass_values_ptr2;
149 boost::shared_ptr<VectorDouble> mass_dots_ptr1;
150 boost::shared_ptr<VectorDouble> mass_dots_ptr2;
153 boost::shared_ptr<ForcesAndSourcesCore>
null;
158 const double T = M_PI / 2.0;
161 double operator() (
const double x)
const {
176 double operator()(
const double x,
const double y,
const double t)
const {
177 double g = cos(
T * x) * cos(
T * y);
188 const double t)
const {
190 double mx = -
T * sin(
T * x) * cos(
T * y);
191 double my = -
T * cos(
T * x) * sin(
T * y);
207 double operator()(
const double x,
const double y,
const double t)
const {
208 double glap = -2.0 * pow(
T, 2) * cos(
T * x) * cos(
T * y);
218 double operator()(
const double x,
const double y,
const double t)
const {
219 double gdot = cos(
T * x) * cos(
T * y);
230 CHKERR m_field.getInterface(simple_interface);
231 CHKERR simple_interface->getOptions();
232 CHKERR simple_interface->loadFile();
237 std::string flux_field) {
239 CHKERR simple_interface->addDomainField(mass_field,
L2,
242 CHKERR simple_interface->addDomainField(flux_field,
HCURL,
245 CHKERR simple_interface->addBoundaryField(flux_field,
HCURL,
250 CHKERR simple_interface->setFieldOrder(mass_field,
order - 1);
251 CHKERR simple_interface->setFieldOrder(flux_field,
order);
252 CHKERR simple_interface->setFieldOrder(
"ERROR", 0);
259 string name = it->getName();
260 const int id = it->getMeshsetId();
261 if (name.compare(0, 14,
"REGION1") == 0) {
266 }
else if (name.compare(0, 14,
"REGION2") == 0) {
271 }
else if (name.compare(0, 14,
"REGION3") == 0) {
276 }
else if (name.compare(0, 14,
"REGION4") == 0) {
281 }
else if (name.compare(0, 14,
"REGION5") == 0) {
293 std::string
internal) {
296 string name = it->getName();
297 if (name.compare(0, 14, natural) == 0) {
299 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
300 natural_bdry_ents,
true);
301 }
else if (name.compare(0, 14, essential) == 0) {
302 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
303 essential_bdry_ents,
true);
304 }
else if (name.compare(0, 14,
internal) == 0) {
305 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
306 internal_edge_ents,
true);
317 block_id,
BLOCKSET, 2, surface,
true);
323 boost::shared_ptr<VectorDouble> &mass_ptr) {
325 vol_ele_slow_rhs->getOpPtrVector().push_back(
331 std::string flux_field,
332 boost::shared_ptr<PreviousData> &data) {
335 vol_ele_slow_rhs->getOpPtrVector().push_back(
338 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
342 internal_edge_ents));
344 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
347 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
353 boost::shared_ptr<PreviousData> &data) {
356 vol_ele->getOpPtrVector().push_back(
360 vol_ele->getOpPtrVector().push_back(
368 boost::shared_ptr<VectorDouble> &mass_ptr,
369 boost::shared_ptr<MatrixDouble> &flux_ptr,
370 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
371 boost::shared_ptr<VectorDouble> &flux_div_ptr) {
375 vol_ele_stiff_rhs->getOpPtrVector().push_back(
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
381 vol_ele_stiff_rhs->getOpPtrVector().push_back(
384 vol_ele_stiff_rhs->getOpPtrVector().push_back(
390 std::string flux_field,
391 boost::shared_ptr<PreviousData> &data,
394 vol_ele_stiff_rhs->getOpPtrVector().push_back(
397 vol_ele_stiff_rhs->getOpPtrVector().push_back(
405 boost::shared_ptr<VectorDouble> &mass_ptr,
406 boost::shared_ptr<MatrixDouble> &flux_ptr) {
408 vol_ele_stiff_lhs->getOpPtrVector().push_back(
411 vol_ele_stiff_lhs->getOpPtrVector().push_back(
417 std::string flux_field,
418 boost::shared_ptr<PreviousData> &data,
421 vol_ele_stiff_lhs->getOpPtrVector().push_back(
424 vol_ele_stiff_lhs->getOpPtrVector().push_back(
427 vol_ele_stiff_lhs->getOpPtrVector().push_back(
430 vol_ele_stiff_lhs->getOpPtrVector().push_back(
437 auto vol_rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
438 vol_ele_slow_rhs->getRuleHook = vol_rule;
439 natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
441 vol_ele_stiff_rhs->getRuleHook = vol_rule;
443 vol_ele_stiff_lhs->getRuleHook = vol_rule;
448 boost::shared_ptr<FaceEle> &initial_ele) {
450 initial_ele->getOpPtrVector().push_back(
458 "SimpleProblem", flux_field, essential_bdry_ents);
498 CHKERR TSSetType(ts, TSARKIMEX);
499 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
502 vol_ele_stiff_lhs,
null,
null);
506 vol_ele_stiff_rhs,
null,
null);
509 vol_ele_slow_rhs,
null,
null);
511 natural_bdry_ele_slow_rhs,
null,
null);
518 std::string flux_field) {
520 post_proc->addFieldValuesPostProc(mass_field);
521 post_proc->addFieldValuesPostProc(flux_field);
528 monitor_ptr,
null,
null);
540 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
541 CHKERR TSSetSolution(ts, X);
542 CHKERR TSSetFromOptions(ts);
546 CHKERR TSGetSNES(ts, &snes);
548 CHKERR SNESGetKSP(snes, &ksp);
550 CHKERR KSPGetPC(ksp, &pc);
551 PetscBool is_pcfs = PETSC_FALSE;
552 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
555 if (is_pcfs == PETSC_TRUE) {
560 problem_ptr->
getName(),
ROW,
"MASS1", 0, 1, &is_mass);
562 problem_ptr->
getName(),
ROW,
"FLUX1", 0, 1, &is_flux);
566 CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
567 CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
570 CHKERR ISDestroy(&is_flux);
571 CHKERR ISDestroy(&is_mass);
585 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
586 CHKERR add_fe(
"MASS1",
"FLUX1");
587 if (nb_species == 2 || nb_species == 3) {
588 CHKERR add_fe(
"MASS2",
"FLUX2");
589 if (nb_species == 3) {
590 CHKERR add_fe(
"MASS3",
"FLUX3");
595 CHKERR simple_interface->setUp();
597 CHKERR set_blockData(material_blocks);
599 CHKERR extract_bd_ents(
"ESSENTIAL",
"NATURAL",
"INTERNAL");
601 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
602 CHKERR extract_initial_ents(2, inner_surface1);
603 CHKERR update_slow_rhs(
"MASS1", mass_values_ptr1);
604 if (nb_species == 1) {
606 "MASS1", data1, data1, data1, material_blocks));
607 }
else if (nb_species == 2 || nb_species == 3) {
608 CHKERR extract_initial_ents(3, inner_surface2);
609 CHKERR update_slow_rhs(
"MASS2", mass_values_ptr2);
610 if (nb_species == 2) {
612 "MASS1", data1, data2, data2, material_blocks));
613 }
else if (nb_species == 3) {
614 CHKERR extract_initial_ents(4, inner_surface3);
615 CHKERR update_slow_rhs(
"MASS3", mass_values_ptr3);
617 "MASS1", data1, data2, data3, material_blocks));
621 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
624 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
625 CHKERR push_slow_rhs(
"MASS1",
"FLUX1", data1);
626 if (nb_species == 2 || nb_species == 3) {
627 CHKERR push_slow_rhs(
"MASS2",
"FLUX2", data2);
628 if (nb_species == 3) {
629 CHKERR push_slow_rhs(
"MASS3",
"FLUX3", data3);
634 CHKERR update_vol_fe(vol_ele_stiff_rhs, data1);
636 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
637 CHKERR update_stiff_rhs(
"MASS1",
"FLUX1", mass_values_ptr1,
638 flux_values_ptr1, mass_dots_ptr1, flux_divs_ptr1);
639 CHKERR push_stiff_rhs(
"MASS1",
"FLUX1", data1,
641 if (nb_species == 2 || nb_species == 3) {
642 CHKERR update_stiff_rhs(
"MASS2",
"FLUX2", mass_values_ptr2,
643 flux_values_ptr2, mass_dots_ptr2, flux_divs_ptr2);
644 CHKERR push_stiff_rhs(
"MASS2",
"FLUX2", data2, material_blocks);
645 if (nb_species == 3) {
646 CHKERR update_stiff_rhs(
"MASS3",
"FLUX3", mass_values_ptr3,
647 flux_values_ptr3, mass_dots_ptr3,
649 CHKERR push_stiff_rhs(
"MASS3",
"FLUX3", data3, material_blocks);
654 CHKERR update_vol_fe(vol_ele_stiff_lhs, data1);
656 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
657 CHKERR update_stiff_lhs(
"MASS1",
"FLUX1", mass_values_ptr1,
659 CHKERR push_stiff_lhs(
"MASS1",
"FLUX1", data1,
661 if (nb_species == 2 || nb_species == 3) {
662 CHKERR update_stiff_lhs(
"MASS2",
"FLUX2", mass_values_ptr2,
664 CHKERR push_stiff_lhs(
"MASS2",
"FLUX2", data2, material_blocks);
665 if (nb_species == 3) {
666 CHKERR update_stiff_lhs(
"MASS3",
"FLUX3", mass_values_ptr3,
668 CHKERR push_stiff_lhs(
"MASS3",
"FLUX3", data3, material_blocks);
676 CHKERR set_integration_rule();
677 dm = simple_interface->getDM();
679 boost::shared_ptr<FaceEle> initial_mass_ele(
new FaceEle(m_field));
681 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
682 CHKERR apply_IC(
"MASS1", inner_surface1,
684 if (nb_species == 2 || nb_species == 3) {
685 CHKERR apply_IC(
"MASS2", inner_surface2, initial_mass_ele);
686 if (nb_species == 3) {
687 CHKERR apply_IC(
"MASS3", inner_surface3, initial_mass_ele);
694 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
696 if (nb_species == 2 || nb_species == 3) {
698 if (nb_species == 3) {
705 post_proc->generateReferenceElementMesh();
709 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
710 CHKERR post_proc_fields(
"MASS1",
"FLUX1");
711 post_proc->addFieldValuesPostProc(
"ERROR");
712 if (nb_species == 2 || nb_species == 3) {
713 CHKERR post_proc_fields(
"MASS2",
"FLUX2");
714 if (nb_species == 3) {
715 CHKERR post_proc_fields(
"MASS3",
"FLUX3");
768 monitor_ptr = boost::shared_ptr<Monitor>(
769 new Monitor(cOmm, rAnk, dm, post_proc, global_error));
775 int main(
int argc,
char *argv[]) {
782 DMType dm_name =
"DMMOFEM";
MoFEM::FaceElementForcesAndSourcesCore FaceEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
int main(int argc, char *argv[])
FaceElementForcesAndSourcesCore BoundaryEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateJacForFaceImpl< 2 > OpCalculateJacForFace
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)
OpSetContravariantPiolaTransformOnEdge OpSetContrariantPiolaTransformOnEdge
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
OpSetContravariantPiolaTransformFaceImpl< 2 > OpSetContravariantPiolaTransformFace
DeprecatedCoreInterface Interface
double operator()(const double x) const
double operator()(const double x, const double y, const double t) const
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
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.
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Calculate divergence of vector field.
Get vector field for H-div approximation.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
keeps basic data about problem
std::string getName() const
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
boost::shared_ptr< VectorDouble > flux_divs_ptr3
MoFEMErrorCode update_stiff_rhs()
MoFEMErrorCode apply_BC(std::string flux_field)
MoFEMErrorCode push_slow_rhs()
boost::shared_ptr< VectorDouble > mass_values_ptr3
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
MoFEMErrorCode post_proc_fields()
MoFEMErrorCode run_analysis()
boost::shared_ptr< MatrixDouble > flux_values_ptr3
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
RDProblem(MoFEM::Core &core, const int order)
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode push_stiff_rhs()
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
boost::shared_ptr< VectorDouble > mass_dots_ptr3
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr)
MoFEMErrorCode setup_system()
MoFEMErrorCode output_result()
MoFEMErrorCode update_stiff_lhs()
boost::shared_ptr< PreviousData > data3
MoFEMErrorCode update_stiff_rhs(std::string mass_field, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr, boost::shared_ptr< VectorDouble > &mass_dot_ptr, boost::shared_ptr< VectorDouble > &flux_div_ptr)
MoFEMErrorCode set_integration_rule()