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[]) {
const string default_options = "-ksp_type gmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-ksp_monitor\n";
string param_file = "param_file.petsc";
if (!static_cast<bool>(ifstream(param_file))) {
std::ofstream file(param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file << default_options;
file.close();
}
}
try {
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 DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm, dm_name);
DataAtIntegrationPts commonData(m_field);
CHKERR commonData.getParameters();
boost::shared_ptr<FEMethod> nullFE;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhs(
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhs(
for (auto &sit : commonData.setOfBlocksData) {
feLhs->getOpPtrVector().push_back(
feLhs->getOpPtrVector().push_back(
feLhs->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
commonData.gradDispPtr));
{{"P", commonData.pPtr}},
{{"U", u_ptr}},
{},
{}));
sit.second));
}
Mat Aij;
{
CHKERR VecGhostUpdateBegin(
d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
d, INSERT_VALUES, SCATTER_FORWARD);
}
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;
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++) {
&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++) {
&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++) {
&fit->second->getLoopFe());
}
}
CHKERR VecAssemblyBegin(F_ext);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetOperators(solver, Aij, Aij);
PetscPrintf(PETSC_COMM_WORLD, "Output file: %s\n", "out.h5m");
CHKERR post_proc.writeFile(
"out.h5m");
}
return 0;
}