#include <MethodForForceScaling.hpp>
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
PetscInt nb_sub_steps = 10;
PetscBool flg_file = PETSC_FALSE;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_atom_test = PETSC_FALSE;
{
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Minimal surface area config", "none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"", 2,
CHKERR PetscOptionsInt(
"-my_nb_sub_steps",
"number of substeps",
"", 10,
&nb_sub_steps, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_partitioned",
"set if mesh is partitioned (this result that "
"each process keeps only part of the mesh",
"", PETSC_FALSE, &is_partitioned, PETSC_NULL);
"-my_is_atom_test",
"is used with testing, exit with error when diverged", "",
PETSC_FALSE, &is_atom_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=BCAST_DELETE;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
} else {
const char *option;
option = "";
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
}
bit_level0.set(0);
0, 2, bit_level0);
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
{
CHKERR moab.get_entities_by_type(root_set, MBTRI, tris,
false);
CHKERR skin.find_skin(root_set, tris,
false, bc_edges);
CHKERR moab.get_connectivity(bc_edges, bc_nodes);
bc_ents.merge(bc_edges);
bc_ents.merge(bc_nodes);
"BC_ELEMENT");
}
bc_ents);
{
DM bc_dm;
CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
CHKERR DMSetType(bc_dm,
"MoFEM");
CHKERR DMSetFromOptions(bc_dm);
CHKERR DMCreateGlobalVector(bc_dm, &T);
minimal_surface_element.
feBcEdge.getOpPtrVector().push_back(
minimal_surface_element.
feBcEdge.getOpPtrVector().push_back(
CHKERR MatAssemblyBegin(
A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
A, MAT_FINAL_ASSEMBLY);
{
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
}
CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
}
DM dm;
{
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm,
"MoFEM");
}
{
CHKERR DMCreateGlobalVector(dm, &T);
}
{
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
"U", mse_common_data, true));
minimal_surface_element.
feRhs.getOpPtrVector().push_back(
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
"U", mse_common_data, false));
minimal_surface_element.
feLhs.getOpPtrVector().push_back(
NULL);
&minimal_surface_element.
feRhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
NULL);
&minimal_surface_element.
feLhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
}
SNESConvergedReason snes_reason;
{
SNES snes;
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR SNESSetFromOptions(snes);
double step_size = 1. / nb_sub_steps;
for (int ss = 1; ss <= nb_sub_steps; ss++) {
CHKERR VecAXPY(T, step_size, T0);
CHKERR SNESSolve(snes, PETSC_NULL, T);
CHKERR SNESGetConvergedReason(snes, &snes_reason);
int its;
CHKERR SNESGetIterationNumber(snes, &its);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"number of Newton iterations = %D\n\n", its);
if (snes_reason < 0) {
if (is_atom_test) {
"atom test diverged");
} else {
break;
}
}
}
CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
}
{
&post_proc);
"PARALLEL=WRITE_PART");
}
{
}
}
return 0;
}