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 <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <PostProcOnRefMesh.hpp>
#include <MinimalSurfaceElement.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();
      CHKERRG(ierr);
 
      // 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
    CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
    // 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
    CHKERR m_field.modify_finite_element_add_field_row(
        "MINIMAL_SURFACE_ELEMENT", "U");
    CHKERR m_field.modify_finite_element_add_field_col(
        "MINIMAL_SURFACE_ELEMENT", "U");
    CHKERR m_field.modify_finite_element_add_field_data(
        "MINIMAL_SURFACE_ELEMENT", "U");
    // Add entities to that element
    CHKERR m_field.add_ents_to_finite_element_by_type(
        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
    CHKERR m_field.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
    MinimalSurfaceElement::CommonData mse_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(
          new MinimalSurfaceElement::OpAssmebleBcRhs("U", F));
      minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
          new MinimalSurfaceElement::OpAssmebleBcLhs("U", A));
      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(
          new MinimalSurfaceElement::OpCalculateCoefficientsAtGaussPts(
              "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(
          new MinimalSurfaceElement::OpCalculateCoefficientsAtGaussPts(
              "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 post_proc.generateReferenceElementMesh();
      CHKERR post_proc.addFieldValuesPostProc("U");
      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);
    }
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}