11static char help[] =
"...\n\n";
32 return 200 * sin(x * 10.) * cos(y * 10.);
38 return sin(x * 10.) * cos(y * 10.);
58 : domainField(
"U"), mField(m_field) {}
98 boundaryMarker = bc_mng->getMergedBlocksMarker({
"BOUNDARY_CONDITION"});
101 bc_mng->getMergedBlocksRange({
"BOUNDARY_CONDITION"});
116 pipeline_mng->getOpDomainLhsPipeline(), {H1});
117 pipeline_mng->getOpDomainLhsPipeline().push_back(
119 pipeline_mng->getOpDomainLhsPipeline().push_back(
121 pipeline_mng->getOpDomainLhsPipeline().push_back(
127 pipeline_mng->getOpDomainRhsPipeline().push_back(
129 pipeline_mng->getOpDomainRhsPipeline().push_back(
131 pipeline_mng->getOpDomainRhsPipeline().push_back(
137 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
139 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
141 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
147 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
149 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
151 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
165 auto domain_rule_lhs = [](int, int,
int p) ->
int {
return 2 * (
p - 1); };
166 auto domain_rule_rhs = [](int, int,
int p) ->
int {
return 2 * (
p - 1); };
168 CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
170 auto boundary_rule_lhs = [](int, int,
int p) ->
int {
return 2 *
p; };
171 auto boundary_rule_rhs = [](int, int,
int p) ->
int {
return 2 *
p; };
172 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
173 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
185 auto ksp_solver = pipeline_mng->
createKSP();
186 CHKERR KSPSetFromOptions(ksp_solver);
196 CHKERR KSPGetPC(ksp_solver, &pc);
197 PetscBool is_pcfs = PETSC_FALSE;
198 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
203 if (is_pcfs == PETSC_TRUE) {
204 IS is_domain, is_boundary;
211 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
212 CHKERR ISDestroy(&is_boundary);
216 CHKERR KSPSetUp(ksp_solver);
222 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
223 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
236 pipeline_mng->getBoundaryLhsFE().reset();
237 pipeline_mng->getDomainRhsFE().reset();
238 pipeline_mng->getBoundaryRhsFE().reset();
240 auto d_ptr = boost::make_shared<VectorDouble>();
241 auto l_ptr = boost::make_shared<VectorDouble>();
245 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(
mField);
247 post_proc_fe->getOpPtrVector().push_back(
249 post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
250 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
251 {{domainField, d_ptr}}, {}, {}, {}));
252 pipeline_mng->getDomainRhsFE() = post_proc_fe;
254 CHKERR pipeline_mng->loopFiniteElements();
255 CHKERR post_proc_fe->writeFile(
"out_result.h5m");
274int main(
int argc,
char *argv[]) {
283 DMType dm_name =
"DMMOFEM";
287 moab::Core mb_instance;
288 moab::Interface &moab = mb_instance;
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#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 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.
auto createDMVector(DM dm)
Get smart vector from DM.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark block DOFs.
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.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
keeps basic data about problem
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 addBoundaryField(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 boundary.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static double boundaryFunction(const double x, const double y, const double z)
Poisson2DNonhomogeneous(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode boundaryCondition()
[Setup problem]
Range boundaryEntitiesForFieldsplit
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode runProgram()
MoFEM::Interface & mField
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode assembleSystem()
[Boundary condition]