Minimal surface area.
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
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);
CHKERR DMRegister_MoFEM(
"MoFEM");
{
DM bc_dm;
CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
CHKERR DMSetType(bc_dm,
"MoFEM");
CHKERR DMMoFEMCreateMoFEM(bc_dm, &m_field,
"DM_BC", bit_level0);
CHKERR DMSetFromOptions(bc_dm);
CHKERR DMMoFEMAddElement(bc_dm,
"BC_ELEMENT");
CHKERR DMCreateGlobalVector(bc_dm, &
T);
minimal_surface_element.
feBcEdge.getOpPtrVector().push_back(
minimal_surface_element.
feBcEdge.getOpPtrVector().push_back(
CHKERR DMoFEMLoopFiniteElements(bc_dm,
"BC_ELEMENT",
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);
CHKERR DMoFEMMeshToGlobalVector(bc_dm,
T, INSERT_VALUES, SCATTER_REVERSE);
}
DM dm;
{
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm,
"MoFEM");
CHKERR DMMoFEMCreateMoFEM(dm, &m_field,
"DM_MINIMAL_AREA", bit_level0);
CHKERR DMMoFEMAddElement(dm,
"MINIMAL_SURFACE_ELEMENT");
}
{
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);
CHKERR DMMoFEMSNESSetFunction(dm,
"MINIMAL_SURFACE_ELEMENT",
&minimal_surface_element.
feRhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
NULL);
CHKERR DMMoFEMSNESSetJacobian(dm,
"MINIMAL_SURFACE_ELEMENT",
&minimal_surface_element.
feLhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
}
SNESConvergedReason snes_reason;
{
SNES snes;
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
CHKERR SNESSetFunction(snes,
F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes,
A,
A, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
Vec T0;
CHKERR DMoFEMMeshToLocalVector(dm, T0, INSERT_VALUES, SCATTER_FORWARD);
double step_size = 1. / nb_sub_steps;
for (int ss = 1; ss <= nb_sub_steps; ss++) {
CHKERR DMoFEMMeshToLocalVector(dm,
T, INSERT_VALUES, SCATTER_REVERSE);
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);
CHKERR DMoFEMMeshToGlobalVector(dm,
T, INSERT_VALUES, SCATTER_REVERSE);
}
{
CHKERR DMoFEMLoopFiniteElements(dm,
"MINIMAL_SURFACE_ELEMENT",
&post_proc);
"PARALLEL=WRITE_PART");
}
{
}
}
return 0;
}
Implementation of minimal surface element.
Post-process fields on refined mesh.
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_ATOM_TEST_INVALID
#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.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Minimal surface equation.
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
Keep date between operators.
Integrate vector on lhs,.
Integrate vector on rhs,.
Evaluate function values and gradients at Gauss Pts.
Evaluate function values and gradients at Gauss Pts.
Implementation of minimal area element.
EdgeElement feBcEdge
Used to calculate dofs on boundary.
SurfaceElement feRhs
To calculate right hand side.
SurfaceElement feLhs
To calculate left hand side.
virtual moab::Interface & get_moab()=0
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode generateReferenceElementMesh()
moab::Interface & postProcMesh