v0.8.15
minimal_surface_area.cpp

Minimal surface area

Implementation is based on: https://www.dealii.org/developer/doxygen/deal.II/step_15.html Look for formulation details on deal.II web page.

Todo:

Make this example with heterogenous approx. apace

Make this example with several refined meshes

/** \file minimal_surface_area.cpp
\example minimal_surface_area.cpp
\brief Minimal surface area
\ingroup minimal_surface_area
Implementation is based on:
<https://www.dealii.org/developer/doxygen/deal.II/step_15.html> Look for
formulation details on deal.II web page.
\todo Make this example with heterogenous approx. apace
\todo Make this example with several refined meshes
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
#include <DirichletBC.hpp>
using namespace MinimalSurfaceEquation;
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
char mesh_file_name[255];
PetscInt order = 2;
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",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
&order, PETSC_NULL);
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);
CHKERR PetscOptionsBool(
"-my_is_atom_test",
"is used with testing, exit with error when diverged", "",
PETSC_FALSE, &is_atom_test, PETSC_NULL);
ierr = PetscOptionsEnd();
// Reade parameters from line command
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE) {
// Read mesh to MOAB
const char *option;
option = "PARALLEL=BCAST_DELETE;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR moab.load_file(mesh_file_name, 0, option);
} else {
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
}
MoFEM::Interface &m_field = core;
// meshset consisting all entities in mesh
EntityHandle root_set = moab.get_root_set();
// Seed all mesh entities to MoFEM database, those entities can be
// potentially used as finite elements or as entities which carry some
// approximation field.
BitRefLevel bit_level0;
bit_level0.set(0);
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 2, bit_level0);
// Add field
// Add entities to field (root_mesh, i.e. on all mesh entities fields are
// approx.) On those entities and lower dimension entities approximation is
// spaned,
CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "U");
// Set approximation oder
CHKERR m_field.set_field_order(root_set, MBTRI, "U", order);
CHKERR m_field.set_field_order(root_set, MBEDGE, "U", order);
CHKERR m_field.set_field_order(root_set, MBVERTEX, "U", 1);
CHKERR m_field.build_fields();
// Add finite element (this defines element declaration comes later)
CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT");
// Set on which fields that element operates
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
// Add entities to that element
root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
CHKERR m_field.add_finite_element("BC_ELEMENT");
CHKERR m_field.modify_finite_element_add_field_row("BC_ELEMENT", "U");
CHKERR m_field.modify_finite_element_add_field_col("BC_ELEMENT", "U");
CHKERR m_field.modify_finite_element_add_field_data("BC_ELEMENT", "U");
// Add entities to that element
Range bc_edges, bc_ents;
{
// Take a skin form the mesh
Skinner skin(&m_field.get_moab());
Range tris;
CHKERR moab.get_entities_by_type(root_set, MBTRI, tris, false);
Range bc_edges;
CHKERR skin.find_skin(root_set, tris, false, bc_edges);
Range bc_nodes;
CHKERR moab.get_connectivity(bc_edges, bc_nodes);
bc_ents.merge(bc_edges);
bc_ents.merge(bc_nodes);
// Add entities to that element
CHKERR m_field.add_ents_to_finite_element_by_type(bc_edges, MBEDGE,
"BC_ELEMENT");
}
// build finite elements
// build adjacencies between elements and degrees of freedom
CHKERR m_field.build_adjacencies(bit_level0);
// Create minimal surface area and elements
MinimalSurfaceElement minimal_surface_element(m_field);
// Data structure for common data
// This class is used to fix dofs on entities
DirichletFixFieldAtEntitiesBc fix_edges_ents(m_field, "U", NULL, NULL, NULL,
bc_ents);
// Solve BC problem
{
// define problems
DM bc_dm;
// Register DM problem
CHKERR DMCreate(PETSC_COMM_WORLD, &bc_dm);
CHKERR DMSetType(bc_dm, "MoFEM");
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(bc_dm, &m_field, "DM_BC", bit_level0);
// Configure DM form line command options (DM itself, solvers,
// pre-conditioners, ... )
CHKERR DMSetFromOptions(bc_dm);
// Add elements to dm (only one here)
CHKERR DMMoFEMAddElement(bc_dm, "BC_ELEMENT");
// Set up problem (number dofs, partition mesh, etc.)
CHKERR DMSetUp(bc_dm);
Mat A;
Vec T, F;
CHKERR DMCreateGlobalVector(bc_dm, &T);
CHKERR VecDuplicate(T, &F);
CHKERR DMCreateMatrix(bc_dm, &A);
// Boundary problem
minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
CHKERR DMoFEMLoopFiniteElements(bc_dm, "BC_ELEMENT",
&minimal_surface_element.feBcEdge);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
{
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetOperators(solver, A, A);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
CHKERR KSPSolve(solver, F, T);
CHKERR KSPDestroy(&solver);
}
// Scatter data on mesh
CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(bc_dm, T, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecDestroy(&T);
CHKERR VecDestroy(&F);
CHKERR MatDestroy(&A);
CHKERR DMDestroy(&bc_dm);
}
// //define problems
// CHKERR m_field.add_problem("DM_MINIMAL_AREA");
// //set refinement level for problem
// CHKERR
// m_field.modify_problem_ref_level_add_bit("DM_MINIMAL_AREA",bit_level0);
// Create discrete manager instance
DM dm;
{
// Register problem
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm, "MoFEM");
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "DM_MINIMAL_AREA", bit_level0);
// Get options fron line command
CHKERR DMSetFromOptions(dm);
// Add elements to dm (only one here)
CHKERR DMMoFEMAddElement(dm, "MINIMAL_SURFACE_ELEMENT");
// Set up data structures
CHKERR DMSetUp(dm);
}
// Create matrices and vectors used for analysis
Vec T, F;
Mat A;
{
CHKERR DMCreateGlobalVector(dm, &T);
CHKERR VecDuplicate(T, &F);
CHKERR DMCreateMatrix(dm, &A);
}
// Adding elements to DMSnes
{
// Right hand side operators (run for each entity in element)
// Calculate inverse of jacobian
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new OpCalculateInvJacForFace(mse_common_data.invJac));
// Apply invJac to shape functions
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(mse_common_data.invJac));
// Below operators are specific to minimal area problem
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
minimal_surface_element.feRhs.getOpPtrVector().push_back(
"U", mse_common_data, true));
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new MinimalSurfaceElement::OpAssembleResidaul("U", mse_common_data));
// Left hand side operators (run for each entity in element)
// Calculate inverse of jacobian
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new OpCalculateInvJacForFace(mse_common_data.invJac));
// Apply invJac to shape functions
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(mse_common_data.invJac));
// Below operators are specific to minimal area problem
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new MinimalSurfaceElement::OpGetDataAtGaussPts("U", mse_common_data));
minimal_surface_element.feLhs.getOpPtrVector().push_back(
"U", mse_common_data, false));
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new MinimalSurfaceElement::OpAssembleTangent("U", mse_common_data));
// Rhs
CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
NULL);
CHKERR DMMoFEMSNESSetFunction(dm, "MINIMAL_SURFACE_ELEMENT",
&minimal_surface_element.feRhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
// Lhs
CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, 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;
// Solve problem
{
SNES snes;
SnesCtx *snes_ctx;
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
// CHKERR SNESSetDM(snes,dm);
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 VecDuplicate(T, &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 VecAXPY(T, step_size, T0);
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) {
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"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 MatView(A,PETSC_VIEWER_DRAW_SELF);
// std::string wait;
// std::cin >> wait;
CHKERR VecDestroy(&T0);
CHKERR SNESDestroy(&snes);
}
{
PostProcFaceOnRefinedMesh post_proc(m_field);
CHKERR DMoFEMLoopFiniteElements(dm, "MINIMAL_SURFACE_ELEMENT",
&post_proc);
// Save data on mesh
CHKERR post_proc.postProcMesh.write_file("out.h5m", "MOAB",
"PARALLEL=WRITE_PART");
}
// Clean and destroy
{
CHKERR DMDestroy(&dm);
CHKERR VecDestroy(&T);
CHKERR VecDestroy(&F);
CHKERR MatDestroy(&A);
}
}
return 0;
}