v0.14.0
Loading...
Searching...
No Matches
minimal_surface_area.cpp

Minimal surface area.

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::Core core(moab);
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);
CHKERR DMRegister_MoFEM("MoFEM");
// 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
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(
new OpCalculateHOJacForFace(jac_ptr));
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
// Apply inv_jac_ptr to shape functions
minimal_surface_element.feRhs.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
// 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 OpCalculateHOJacForFace(jac_ptr));
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
// Apply inv_jac_ptr to shape functions
minimal_surface_element.feLhs.getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
// 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);
CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL, 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);
CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, 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;
}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
Implementation of minimal surface element.
Post-process fields on refined mesh.
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
@ F
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.
char mesh_file_name[255]
const double T
Minimal surface equation.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
constexpr AssemblyType A
Fix dofs on entities.
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.
Managing BitRefLevels.
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.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.
MoFEMErrorCode generateReferenceElementMesh()
moab::Interface & postProcMesh