v0.13.2
Loading...
Searching...
No Matches
poisson_3d_homogeneous.cpp

Poisson problem 3D.

Poisson problem 3D

/**
* @file poisson_3d_homogeneous.cpp
* @example poisson_3d_homogeneous.cpp
* @brief Poisson problem 3D
*
* @copyright Copyright (c) 2023
*
*/
using namespace MoFEM;
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
// MoFEM interfaces
// Field name and approximation order
std::string domainField;
int oRder;
};
: domainField("U"), mField(m_field) {}
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
auto bc_mng = mField.getInterface<BcManager>();
// Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
// can use regular expression to put list of blocksets;
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
domainField, true);
// Remove BCs from cubit TEMPERATURESET, i.e. set by cubit, and meshsets named
// FIX_SCALAR (default name to name boundary conditions for scalar fields)
CHKERR bc_mng->removeBlockDOFsOnEntities<TemperatureCubitBcData>(
simpleInterface->getProblemName(), domainField, true, false, true);
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
{ // Push operators to the Pipeline that is responsible for calculating LHS
CHKERR AddHOOps<3, 3, 3>::add(pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating LHS
constexpr int space_dim = 3;
auto set_values_to_bc_dofs = [&](auto &fe) {
auto get_bc_hook = [&]() {
return hook;
};
fe->preProcessHook = get_bc_hook();
};
auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
using DomainEle =
auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
grad_u_vals_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInternal(domainField, grad_u_vals_ptr,
[](double, double, double) constexpr { return -1; }));
};
pipeline_mng->getOpDomainRhsPipeline(), {H1});
set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
calculate_residual_from_set_values_on_bc(
pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
}
auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
auto rule_rhs = [](int, int, int p) -> int { return p; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR KSPSetUp(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = smartCreateDMVector(dm);
// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcVolEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{domainField, u_ptr}},
{},
{},
{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_result.h5m");
}
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
Poisson3DHomogeneous poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
static Index< 'p', 3 > p
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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:511
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcVolEle
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
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.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
Definition of the temperature bc data structure.
Definition: BCData.hpp:301
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode setIntegrationRules()
Poisson3DHomogeneous(MoFEM::Interface &m_field)
MoFEMErrorCode boundaryCondition()