Main implementation of U-P (mixed) finite element.
Main implementation of U-P (mixed) finite element.
using namespace boost::numeric;
static char help[] =
"-my_order_p approximation order_p \n"
"-my_order_u approximation order_u \n"
"-is_partitioned is_partitioned? \n";
int main(
int argc,
char *argv[]) {
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-ksp_monitor\n";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
auto moab_comm_wrap =
boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
PetscBool flg_file;
int order_p = 2;
int order_u = 3;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool flg_test = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Mix elastic problem",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-my_order_p",
"approximation order_p",
"", order_p,
&order_p, PETSC_NULL);
CHKERR PetscOptionsInt(
"-my_order_u",
"approximation order_u",
"", order_u,
&order_u, PETSC_NULL);
CHKERR PetscOptionsBool(
"-is_partitioned",
"is_partitioned?",
"",
is_partitioned, &is_partitioned, PETSC_NULL);
CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE) {
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
} else {
const char *option;
option = "";
}
bit_level0.set(0);
0, 3, bit_level0);
3);
{
"MESH_NODE_POSITIONS");
}
"MESH_NODE_POSITIONS");
DM dm;
DMType dm_name = "DM_ELASTIC_MIX";
CHKERR DMRegister_MoFEM(dm_name);
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm, dm_name);
CHKERR DMMoFEMCreateMoFEM(dm, &m_field, dm_name, bit_level0);
CHKERR DMMoFEMAddElement(dm,
"ELASTIC");
CHKERR DMMoFEMAddElement(dm,
"PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm,
"FORCE_FE");
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
boost::shared_ptr<FEMethod> nullFE;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhs(
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhs(
feLhs->getOpPtrVector().push_back(
feLhs->getOpPtrVector().push_back(
feLhs->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc.getOpPtrVector().push_back(
post_proc.getOpPtrVector().push_back(
post_proc.getOpPtrVector().push_back(
post_proc.getOpPtrVector().push_back(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
{{"P", commonData.pPtr}},
{{"U", u_ptr}},
{},
{}));
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), commonData,
sit.second));
}
Mat Aij;
Vec d, F_ext;
{
CHKERR DMCreateGlobalVector_MoFEM(dm, &d);
CHKERR VecDuplicate(d, &F_ext);
CHKERR DMCreateMatrix_MoFEM(dm, &Aij);
CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
}
feLhs->ksp_B = Aij;
feLhs->ksp_f = F_ext;
boost::shared_ptr<DirichletDisplacementBc> dirichlet_bc_ptr(
dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
CHKERR DMoFEMLoopFiniteElements(dm,
"ELASTIC", feLhs);
boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
F_ext, "U");
{
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumann_forces.begin();
for (; mit != neumann_forces.end(); mit++) {
CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
&mit->second->getLoopFe());
}
}
boost::ptr_map<std::string, NodalForce> nodal_forces;
{
boost::ptr_map<std::string, NodalForce>::iterator fit =
nodal_forces.begin();
for (; fit != nodal_forces.end(); fit++) {
CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
&fit->second->getLoopFe());
}
}
boost::ptr_map<std::string, EdgeForce> edge_forces;
{
boost::ptr_map<std::string, EdgeForce>::iterator fit =
edge_forces.begin();
for (; fit != edge_forces.end(); fit++) {
CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
&fit->second->getLoopFe());
}
}
CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
CHKERR DMMoFEMKSPSetComputeOperators(dm,
"ELASTIC", feLhs, nullFE, nullFE);
CHKERR VecAssemblyBegin(F_ext);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetOperators(solver, Aij, Aij);
CHKERR KSPSolve(solver, F_ext, d);
CHKERR DMoFEMMeshToGlobalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
CHKERR DMoFEMLoopFiniteElements(dm,
"ELASTIC", &post_proc);
PetscPrintf(PETSC_COMM_WORLD, "Output file: %s\n", "out.h5m");
CHKERR post_proc.writeFile(
"out.h5m");
}
return 0;
}
const std::string default_options
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
boost::shared_ptr< MatrixDouble > gradDispPtr
std::map< int, BlockData > setOfBlocksData
boost::shared_ptr< VectorDouble > pPtr
MoFEMErrorCode getParameters()
Set Dirichlet boundary conditions on displacements.
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Get value at integration points for scalar field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.