Soap film spanned on wire

Table of Contents

• Theory
• Installation
  • Updating module
• Definition and Declarations
  • Declarations & definition
    • Declaring and definition of field
    • Declaring & definition finite elements
  • Declaring and defining problem and registering discrete manager
  • Solving problems
    • Solving linear boundary problem
    • Solving nonlinear minimal area equation
• Implementing finite elements
  • Finite element on boundary
  • Finite element for minimal surface area problem
• Example
• Future work and contribution
• The plain program

Example of implementation, step by step tutorial.

This documents assumes that reader is familiar with hierarchical approximation spaces, if not pleas look into paper [1]. In such approximation degrees of freedom are not only associated with vertices but edges and other two- and three- dimensional entities. MoFEM exploit that property and use it as a tool to make efficient and well structured code.

Soap film on rectangular wire loop with g(x,y)=sin(Pi*x)

Theory

This problem formulation follows tutorial from deal.II https://www.dealii.org/developer/doxygen/deal.II/step_15.html. This problem describe the surface spanned by soap film enclosed between a rigid wire. For simplicity shape of a wire is given by curve \(z=g(x,y)\).

Problem equation in the strong form is given by

\[ -\nabla \cdot \left( \frac{1}{\sqrt{1+\|u\|^2}}\nabla u \right) = 0 \; \textrm{in}\,\Omega,\\ u = g(x,y) \; \textrm{on}\,\partial\Omega \]

Above equation is linearized and applying standard procedure of integration by parts expressed finally in the week form

\[ \int_\Omega \nabla \phi \cdot \left\{ a_i \nabla \delta u_{i+1} - a_i^3 \nabla u_i \cdot \nabla \delta u_{i+1} \nabla u_i \right\} \textrm{d}\Omega + \int_\Omega \nabla \phi \cdot \left\{ a_i \nabla u_{i+1} \right\} \textrm{d}\Omega = 0 \]

where

\[ a_i = \frac{1}{\sqrt{1+\|u_i\|^2}} \]

and \(i\) is iteration number where function values at step \(s\) are updated us follows

\[ u^{s+1} = u^s + \Delta u_{i+1},\\ \Delta u_{i+1} = \Delta u_i + \delta u_{i+1}. \]

Installation

MoFEM module (like MatLab toolkit) is independent implementation, which own repository, different owners and license. Adding module is straightforward, CMake search in subdirectories for a file InstalledAddModule.cmake, directory in which such file is placed become an user module. You can look to the How to add user module? for additional information about module installation and its directory structure.

It is good practice that each module has tests verifying correctness of the installation. Having tests is essential for code development and ensures compatibility of a module which MoFEM library and other modules.

The minimal surface problem is implemented as a stand alone independent module which own repository and GNU Lesser General Public License. A module is installed from Bitbucket repository (https://bitbucket.org/likask/mofem_um_minimal_surface_equation) by command,

export MOFEM_INSTALL_DIR=$HOME/mofem_installation
git clone https://bitbucket.org/likask/mofem_um_minimal_surface_equation \
$MOFEM_INSTALL_DIR/mofem-cephas/mofem/users_modules/minimal_surface_equation

Modules is compiled with commands

cd $MOFEM_INSTALL_DIR/users_modules
touch CMakeCache.txt && make

Finally installation is tested with

ctest -R minimal_surface_equation_test -VV

Updating module

Since modules are located in different independent repositories, for connivence update of all installed modules is automatized. You can do that with line command;

cd $MOFEM_INSTALL_DIR/users_modules
make update_users_modules && make

you can as well update individual module, for example

cd $MOFEM_INSTALL_DIR/users_modules
make update_minimal_surface_equation && make

Alternatively to update this module you can simply do;

cd $MOFEM_INSTALL_DIR/users_modules/minimal_surface_equation
git pull && make

Definition and Declarations

In MoFEM like in programming languages like C++ we use declarations, definitions and implementation. A definition of finite element is what MoFEM needs to build data structures in order to make vectors, matrices and set up solvers. A implementation of finite element is what MoFEM needs to do calculations, i.e. calculate element matrices and vectors. Such approach allow to have the same declaration but several definitions (implementations) of the same element.

Declarations & definition

Declaring and definition of field

In this problem we have only one scalar field of a soap surface displacements in \(z\) direction. To solve problem using week form shown above we demand that approximation field belongs to H1 space (i.e. conforming approximation). Since problem is 2D, we span that field on triangles, edges and nodes;

// 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(); 

Note that initially we set the same approximation order to each entity. However in general approximation order can be heterogenous. Approximation order of vertices is always "1" for H1 space, since for given nodes values uniquely piecewise linear function is given.

If it is needed additional abstraction, the field can carry additional information about its covariant and contravariant base, see MoFEM::CoordSys. We do not need this here, since we work only with one component of displacement field in Cartesian base.

Declaring & definition finite elements

In MoFEM finite element is entity (Vertex, Edge, Triangle, Quad, Tet, Prism ... Meshset), on which we do some calculations, f.e. integrate internal force vector.

Definition of a finite element resembles how we write bilinear form, f.e. problem in the week form is given by general equation

\[ b(v,u;\alpha) = f(v)\;\in\Omega \]

where \(u,v\in\mathcal{V} \subset H_1(\Omega)\) belongs to some finite vector space and \(\alpha = \alpha(x,y) \in L_2(\Omega)\) is some field of coefficients. A testing function (wighting) space, tested (approximated) field and data (coefficient) field on which bilinear form operates, need to be specified when finite element is declared. Note that in general bilinear form can be decomposed on several parts (by introducing several elements) to split complex problem on simpler tasks.

In this example we declaring two elements, one Edge elements used to calculate approximation of function \(u=g(x,y)\) on boundary \(\partial \Omega\), and another on Triangles to formulate main problem of minimal surface equation in \(\Omega\).

The element in film domain is declared as follows

// Add (declare) finite element (this defines element declaration comes later)
CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT"); 
// Set (define) 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"); 

and for boundary elements element is declared in analogically, as follows

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"); 

In addition we need to set domain of the finite element, it is done by adding entities to the finite element, for MINIMAL_SURFACE_ELEMENT we will use all triangle elements in the root meshset;

// Add entities to that element
CHKERR m_field.add_ents_to_finite_element_by_type(root_set,MBTRI,"MINIMAL_SURFACE_ELEMENT"); 

For BC_ELEMENT, we will take a skin form triangles to obtain edges on the boundary and add them to the boundary finite element;

// Take a skin form the mesh
Skinner skin(&m_field.get_moab());
Range tris;
rval = moab.get_entities_by_type(root_set,MBTRI,tris,false); 
Range bc_edges;
rval = skin.find_skin(root_set,tris,false,bc_edges); 
// Add entities to that element
CHKERR m_field.add_ents_to_finite_element_by_EDGEs(
  bc_edges,"BC_ELEMENT"
); 

Note that dimension of a domain is not explicitly given, it is associated with the type of an entity, in this case we use triangles and edges. However user can do more exotic combinations, mixing different dimensions and do generic element from a meshset. Dimension of a field is implicitly defined not by a finite element, but fields on which it operates and structure MoFEM::CoordSys carried with each field.

Declaring and defining problem and registering discrete manager

The KSP solver, SNES solver or TS solver from PETSc library could be run directly using native MoFEM interface, optionally PETSc Discrete Manager (DM) http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/ can be used by MoFEM DM interface, see Distributed mesh manager.

In this example we will solve sequentially two problems, linear problem finding function approximation on the boundary and nonlinear problem on the film domain. First problem is defined as follows;

DM bc_dm;
DMType dm_name = "DM_BC";
// Register DM problem
CHKERR DMRegister_MoFEM(dm_name); 
CHKERR DMCreate(PETSC_COMM_WORLD,&bc_dm);
CHKERR DMSetType(bc_dm,dm_name);
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(bc_dm,&m_field,dm_name,bit_level0); 
// Configure DM from 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); 

where bit_level0 declares on which refinement level problem is solved. In principle user can have several problems on different mesh refinements, or problem which spanning on several nested meshes.

Analogically second problem is created as follows;

DM dm;
DMType dm_name = "DM_MINIMAL_AREA";
// Register problem
CHKERR DMRegister_MoFEM(dm_name); 
CHKERR DMCreate(PETSC_COMM_WORLD,&dm);
CHKERR DMSetType(dm,dm_name);
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(dm,&m_field,dm_name,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); 

Once DMs are created for each problem, using PETSc interface http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/ we can create vectors and matrices for each of the problems, f.e.

Vec T,F;
Mat A;
// Create vectors
CHKERR DMCreateGlobalVector(dm,&T); 
CHKERR VecDuplicate(T,&F); 
// Create matrix
CHKERR DMCreateMatrix(dm,&A); 

Such implementation make transparent for user how vectors and matrices are created. Since DM is defined for MoFEM problem, problem is composed from set of finite elements and those are using fields, for given space of those fields (L2,H1,Hdiv,Hcurl) and given entities caring those fields, adjacency matrix is uniquely defined.

In MoFEM solution is build "bottom to up", such that fields, approximations, finite elements and their implementation is a problem independent. This enables to use finite element to solve problems, independently on their dimension and number of fields on which those problems operates.

In this example, we exploit this to some extend. The 1d problem is solved first on the edges and then the 2d problem is solved on triangles. Those two problems operates on the same field and since field is independent from problem, data transfer from one to another is effortless.

Solving problems

Solving linear boundary problem

To solve problem the right hand side vector and the left hand side matrix need to be calculated,

// Loop, integrate and assemble
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); 

The element implementation is hidden in its object instance minimal_surface_element.feBcEdge, implementation of this class is discussed in following sections. Having vector and matrices assembled system of equations is solved using KSP solver

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); 

Once problem is solved, results are stored in MoFEM database, independently from particular problem numeration, partitioning on several processors or refinement level

CHKERR DMoFEMMeshToGlobalVector(bc_dm,T,INSERT_VALUES,SCATTER_REVERSE); 

Solving nonlinear minimal area equation

Solution of a nonlinear problem in MoFEM is as simple as solving linear problem, taking advantage from PETSc functionality. To run SNES solver one need to register finite element object instances responsible for calculation of a tangent matrix and a vector of residuals, this is simply done as follows;

//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); 

where fix_edges_ents, minimal_surface_element.feRhs and minimal_surface_element.feLhs are particular finite element implementations. Implementation of fix_edges_ents is not discussed here, it is general propose finite element to enforce essential boundary conditions, in this case fixes DOFs on domain boundary. Implementation of minimal_surface_element.feRhs and minimal_surface_element.feLhs is discussed in following sections.

User can add several elements, for each Newton iteration, elements instances are executed in order in which were added (registered) in MoFEM SNES solver. Low level element class is here MoFEM::FEMethod. Usually user not uses directly this class, it implement using set of so called data entity operators as is shown bellow. Such structure allow to break down complex problem on smaller tasks, easy to debug, reuse in different modules and tests and last but not least produce efficient code.

To create nonlinear solver instance, linked to DM, follow code;

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); 

Note that type of SNES solver, line searchers, parameters are set from line command, see for details https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/, in particular https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESSetFromOptions.html

Analyzed problem in this example strongly nonlinear it is convenient to apply sub-stepping, as in the following code;

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 VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); 
CHKERR VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); 
CHKERR DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); 

where DMoFEMMeshToLocalVector get values from the field (previously calculated in boundary problem) and stores values in vector T0. This vector is scaled by size step size and added to unknown vector.

Implementing finite elements

A finite element is simply understand here as a object doing some calculations (usually integration) on some domain associated with element entity, f.e. example for edge is a segment emerged in 1d/2d/3d space depending on needed level of abstraction, for triangle is a 2d domain (potential emerged in some higher order dimension space). Low level implementation starts with MoFEM::FEMethod, however here we will implement elements with some "driver" classes which will do some typical work for a user. This section as well introduces concept user data operator (see MoFEM::DataOperator). In this implementation two types of elements are used Edge Element and Face Element.

Finite element operates on one or several fields, where approximation spaces of fields are span on some set of entities, specific for each field space and approximation order. For example, space H1 and order 1, span only on vertices, order 2 span on vertices and edges, order 3 span additionally on faces, and so on. Another example, space Hdiv in 3d, span only on faces and volumes, depending on approximation order.

Managing complexities for hierarchical and heterogenous approximation fields on unstructured meshes could be difficult. MoFEM try to make that easy but yet give maximal flexibility. Inside of each element, MoFEM loops over all entities or all possible combinations, according to canonical element numeration [tautges_moab:2004]. For each looped entity ordered set of operators (implemented by user) is executed, where from operator user get easy access to finite element data and data on specific entity, like shape functions, vales of degrees of freedom, indices and other data more specific to finite element entity type.

Implementation of this problem is in MinimalSurfaceElement.hpp particular in class MinimalSurfaceEquation::MinimalSurfaceElement.

Finite element on boundary

Problem boundary is implemented with the MoFEM::EdgeElementForcesAndSourcesCore, see

struct MinimalSurfaceEquation::MinimalSurfaceElement::EdgeElement: public EdgeElementForcesAndSourcesCore {
  EdgeElement(MoFEM::Interface &m_field):
  EdgeElementForcesAndSourcesCore(m_field) {}
  /** \brief Set integrate rule
  */
  int getRule(int order) { return 2*order; };
};
EdgeElement feBcEdge;

where order of appropriate oder of quadrature is set to integrate least square problem

\[ \int_L \mathbf{N}^\textrm{T}\mathbf{N} \textrm{d}L \; \mathbf{u} = \int_L \mathbf{N}^\textrm{T} g(x,y) \textrm{d}L \]

where \(\mathbf{u}\) represents vector of unknown DOFs on boundary of the domain.

The right hand side vector is calculated using instance of

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpAssmebleBcRhs: public MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator {
  Vec vF;
  OpAssmebleBcRhs(const string field_name,Vec v_f):
  EdgeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
  vF(v_f) {}
  virtual double evalFunction(double x,double y) {
    return sin(2*M_PI*(x+y));
  }
  VectorDouble nF;
  PetscErrorCode doWork(
    int side,EntityType type,DataForcesAndSourcesCore::EntData &data
  ) {
    MoFEMFunctionBegin;
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    nF.resize(nb_dofs,false);
    nF.clear();
    int nb_gauss_pts = data.getN().size1();
    for(int gg = 0;gg!=nb_gauss_pts;gg++) {
      double val = getLength()*getGaussPts()(1,gg);
      double x = getCoordsAtGaussPts()(gg,0);
      double y = getCoordsAtGaussPts()(gg,1);
      noalias(nF) += val*evalFunction(x,y)*data.getN(gg);
    }
    CHKERR VecSetValues(
      vF,data.getIndices().size(),&data.getIndices()[0],&nF[0],ADD_VALUES
    ); 
    MoFEMFunctionReturn(0);
  }
};

This UserDataOperator is has type OPROW, it means that it operates on fields which are on rows, "left side" of bilinear form. In that class user has access to range of information on different abstraction levels. Importantly here information about numerical integration, finite element edge length or global position of integration points is accessed. Member function doWork is executed of each entity adjacent to finite element, i.e. vertices and edge. If will be need this class could be inherited and virtual function evalFunction overload, to change function on the boundary of the domain.

The left hand side operator is used to calculate stiffness matrix,

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpAssmebleBcLhs: public MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator {
  Mat mA;
  OpAssmebleBcLhs(const string field_name,Mat m_a):
  EdgeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROWCOL),
  mA(m_a) {}
  MatrixDouble K,transK;
  PetscErrorCode doWork(
    int row_side,int col_side,
    EntityType row_type,EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data
  ) {
    MoFEMFunctionBegin;
 
    int row_nb_dofs = row_data.getIndices().size();
    if(row_nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    int col_nb_dofs = col_data.getIndices().size();
    if(col_nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    K.resize(row_nb_dofs,col_nb_dofs,false);
    K.clear();
    int nb_gauss_pts = row_data.getN().size1();
    for(int gg = 0;gg!=nb_gauss_pts;gg++) {
      double val = getLength()*getGaussPts()(1,gg);
      VectorAdaptor row_n = row_data.getN(gg);
      VectorAdaptor col_n = col_data.getN(gg);
      noalias(K) += val*outer_prod(row_n,col_n);
    }
    CHKERR MatSetValues(
      mA,
      row_nb_dofs,&*row_data.getIndices().data().begin(),
      col_nb_dofs,&*col_data.getIndices().data().begin(),
      &*K.data().begin()
      ,ADD_VALUES
    ); 
    // Matrix is symmetric, assemble off diagonal part
    if(row_type != col_type || row_side != col_side) {
      transK.resize(col_nb_dofs,row_nb_dofs,false);
      noalias(transK) = trans(K);
      CHKERR MatSetValues(
        mA,
        col_nb_dofs,&*col_data.getIndices().data().begin(),
        row_nb_dofs,&*row_data.getIndices().data().begin(),
        &*transK.data().begin()
        ,ADD_VALUES
      ); 
    }
    MoFEMFunctionReturn(0);
  }
};

This operator has type OPROWCOL, it means that it loops over rows and columns and all possible combinations, taking in account that problem is symmetric. The values of shape functions at integration points are obtained by

VectorAdaptor row_n = row_data.getN(gg);
VectorAdaptor col_n = col_data.getN(gg);

where indices of degrees of freedom can be get by

VectorInt &row_indices = row_data.getIndices();
VectorInt &col_indices = col_data.getIndices();

Note that each operator has overload function doWork, depending which variant of a function is overloaded work is done on rows, columns or combinations or possible rows and columns. This member function is run for all adjacent entities to entity of the finite element.

Each finite element keeps vector of smart pointers to user data operators. The vector of smart pointer to operators in this particulate is accessed by

// Boundary problem
MinimalSurfaceElement minimal_surface_element(m_field);
boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &list_of_operators = minimal_surface_element.feBcEdge.getOpPtrVector();

Finally element is composed by adding to it operators,

// Push operators to finite element
list_of_operators.push_back(new MinimalSurfaceElement::OpAssmebleBcRhs("U",F));
list_of_operators.push_back(new MinimalSurfaceElement::OpAssmebleBcLhs("U",A));

which can be executed by Discrete Manager;

CHKERR DMoFEMLoopFiniteElements(bc_dm,"BC_ELEMENT",&minimal_surface_element.feBcEdge); 

Finite element for minimal surface area problem

The minimal surface area is decomposed on several smaller operations. Note that implementation of operator and finite element is independent of the problem and if there will be need, all showed here operators can be reused in different context in to solve different problem.

Implementation starts with data structure (see MinimalSurfaceEquation::MinimalSurfaceElement::CommonData), shared between finite elements and user data operators. It could be understand as small memory pool, used to sharing and data exchange. Is used as well to accumulate results from different entities;

struct MinimalSurfaceEquation::MinimalSurfaceElement::CommonData {
  MatrixDouble gradU; ///< nb_gauss_pts x 2 gradients at integration pts
  VectorDouble normGradU2; ///< size of nb_gauss_pts, norm of gradient
  VectorDouble aN; ///< size of nb_gauss_pts, \f$a_n = \frac{1}{\sqrt{1+\|u_i\|^2}}\f$
  MatrixDouble aNbyGradU; ///< nb_gauss_pts x 2, \f$a_n \nabla u\f$
  MatrixDouble aNpow3byGradU; ///< nb_gauss_pts x 2, \f$a_n^3 \nabla u\f$
};

Next we implement a finite element on face entity, in this particular case on triangle

struct MinimalSurfaceEquation::MinimalSurfaceElement::SurfaceElement: public MoFEM::FaceElementForcesAndSourcesCore {
  SurfaceElement(MoFEM::Interface &m_field):
  FaceElementForcesAndSourcesCore(m_field) {}
  /** \brief Set integrate rule
  */
  int getRule(int order) { return 2*(order-1)+1; };
};

where here only function controlling quadrature rule is overloaded.

The actual work, i.e. integration of residual and tangent matrix, is made by operators. Here we break down work on smaller four task. In first user data operator is implemented to calculate displacement gradient at each integration points;

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpGetDataAtGaussPts: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
  CommonData &commonData;
  OpGetDataAtGaussPts(const string field_name,CommonData &common_data):
  FaceElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPCOL),
  commonData(common_data) {
  }
  PetscErrorCode doWork(
    int side,EntityType type,DataForcesAndSourcesCore::EntData &data
  ) {
    MoFEMFunctionBegin;
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    int nb_gauss_pts = data.getN().size1();
    // Element loops over entities, it start from vertices
    if(type == MBVERTEX) {
      // Setting size of common data vectors
      commonData.gradU.resize(nb_gauss_pts,2);
      commonData.gradU.clear();
    }
    for(int gg = 0;gg!=nb_gauss_pts;gg++) {
      MatrixAdaptor grad_shape_fun = data.getDiffN(gg);
      ublas::matrix_row<MatrixDouble > grad_at_gauss_pt(commonData.gradU,gg);
      noalias(grad_at_gauss_pt) += prod(trans(grad_shape_fun),data.getFieldData());
    }
    MoFEMFunctionReturn(0);
  }
};

This is operator of type OPCOL, it is executed on adjacent entities to finite element entity (in this case, triangle element is adjacent to three vertices, three edges and one face). For each finite element loop over entities starts always from lowest possible entity dimension i.e. vertices, this is exploited here to setup common data structure, i.e. to set sizes and clear vectors and matrixes.

Next coefficients in linearized week from are evaluated;

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpCalculateCoefficientsAtGaussPts: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
  CommonData &commonData;
  bool rhsAndNotLhs;
  OpCalculateCoefficientsAtGaussPts(const string field_name,CommonData &common_data,bool rhs_and_not_lhs):
  FaceElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
  commonData(common_data),
  rhsAndNotLhs(rhs_and_not_lhs) {
  }
  PetscErrorCode doWork(
    int side,EntityType type,DataForcesAndSourcesCore::EntData &data
  ) {
    MoFEMFunctionBegin;
    // Element loops over entities, it start from vertices
    if(type == MBVERTEX) {
      int nb_gauss_pts = data.getN().size1();
      commonData.normGradU2.resize(nb_gauss_pts,false);
      for(int gg = 0;gg!=nb_gauss_pts;gg++) {
        ublas::matrix_row<MatrixDouble > grad_at_gauss_pt(commonData.gradU,gg);
        commonData.normGradU2[gg] = inner_prod(grad_at_gauss_pt,grad_at_gauss_pt);
      }
      commonData.aN.resize(nb_gauss_pts,false);
      for(int gg = 0;gg!=nb_gauss_pts;gg++) {
        commonData.aN[gg] = 1./sqrt(1+commonData.normGradU2[gg]);
      }
      if(rhsAndNotLhs) {
        // needed at right hand side when residual is calculated
        commonData.aNbyGradU.resize(nb_gauss_pts,2,false);
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ublas::matrix_row<MatrixDouble > grad_at_gauss_pt(commonData.gradU,gg);
          for(int rr = 0;rr!=2;rr++) {
            commonData.aNbyGradU(gg,rr) = commonData.aN[gg]*grad_at_gauss_pt[rr];
          }
        }
      } else {
        // needed at left hand side when matrix is calculated
        commonData.aNpow3byGradU.resize(nb_gauss_pts,2,false);
        for(int gg = 0;gg!=nb_gauss_pts;gg++) {
          ublas::matrix_row<MatrixDouble > grad_at_gauss_pt(commonData.gradU,gg);
          for(int rr = 0;rr!=2;rr++) {
            commonData.aNpow3byGradU(gg,rr) = pow(commonData.aN[gg],3)*grad_at_gauss_pt[rr];
          }
        }
      }
    }
    MoFEMFunctionReturn(0);
  }
};

This operator has two variants controlled by variable rhsAndNotLhs, depending if is executed when left or right hand side is calculated, different set of coefficients is needed. Note that coefficients are calculated for given value of gradient already calculated in different data operator.

Once coefficients are calculated residual vector is assembled;

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpAssembleResidaul: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
  CommonData &commonData;
  OpAssembleResidaul(const string field_name,CommonData &common_data):
  FaceElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
  commonData(common_data) {
  }
  VectorDouble rEsidual;
  PetscErrorCode doWork(
    int side,EntityType type,DataForcesAndSourcesCore::EntData &data
  ) {
    MoFEMFunctionBegin;
 
    int nb_dofs = data.getIndices().size();
    if(nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    rEsidual.resize(nb_dofs,false);
    rEsidual.clear();
    int nb_gauss_pts = data.getN().size1();
    for(int gg = 0;gg!=nb_gauss_pts;gg++) {
      double val = getGaussPts()(2,gg)*getArea();
      ublas::matrix_row<MatrixDouble > an_by_grad_u(commonData.aNbyGradU,gg);
      noalias(rEsidual)+=val*prod(data.getDiffN(gg),an_by_grad_u);
    }
    CHKERR VecSetValues(
      getFEMethod()->snes_f,
      data.getIndices().size(),
      &*data.getIndices().data().begin(),
      &rEsidual[0],
      ADD_VALUES
    ); 
    MoFEMFunctionReturn(0);
  }
};

and tangent matrix

struct MinimalSurfaceEquation::MinimalSurfaceElement::OpAssembleTangent: public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
  CommonData &commonData;
  OpAssembleTangent(const string field_name,CommonData &common_data):
  FaceElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROWCOL),
  commonData(common_data) {
    sYmm = false; // matrix is not symmetric
  }
  MatrixDouble kTangent,auxMat;
  VectorDouble auxVec;
  PetscErrorCode    doWork(
    int row_side,int col_side,
    EntityType row_type, EntityType col_type,
    DataForcesAndSourcesCore::EntData &row_data,
    DataForcesAndSourcesCore::EntData &col_data
  ) {
    MoFEMFunctionBegin;
 
    int row_nb_dofs = row_data.getIndices().size();
    if(row_nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    int col_nb_dofs = col_data.getIndices().size();
    if(col_nb_dofs == 0) {
      MoFEMFunctionReturnHot(0);
    }
    kTangent.resize(row_nb_dofs,col_nb_dofs,false);
    kTangent.clear();
    auxVec.resize(col_nb_dofs,false);
    auxMat.resize(col_nb_dofs,2,false);
    int nb_gauss_pts = row_data.getN().size1();
    for(int gg = 0;gg!=nb_gauss_pts;gg++) {
      double val = getGaussPts()(2,gg)*getArea();
      MatrixAdaptor row_diff_n = row_data.getDiffN(gg);
      MatrixAdaptor col_diff_n = col_data.getDiffN(gg);
      double an = commonData.aN[gg];
      ublas::matrix_row<MatrixDouble > grad_u(commonData.gradU,gg);
      ublas::matrix_row<MatrixDouble > an_pow3_by_grad_u(commonData.aNpow3byGradU,gg);
      noalias(auxVec) = prod(col_diff_n,an_pow3_by_grad_u);
      noalias(auxMat) = outer_prod(auxVec,grad_u);
      noalias(kTangent) += val*prod(row_diff_n,trans(an*col_diff_n-auxMat));
    }
    CHKERR MatSetValues(
      getFEMethod()->snes_B,
      row_nb_dofs,&*row_data.getIndices().data().begin(),
      col_nb_dofs,&*col_data.getIndices().data().begin(),
      &*kTangent.data().begin()
      ,ADD_VALUES
    ); 
    MoFEMFunctionReturn(0);
  }
};

Finally when finite element and entity operators are implemented, we add operators to finite elements instances, we do that staring from the right hand side operator

// Data structure for common data
MinimalSurfaceElement::CommonData mse_common_data;
// Right hand side operators (run for each entity in element)
// Calculate inverse of the Jacobian
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 MoFEM::OpSetInvJacH1ForFace(inv_jac_ptr));
// Below operators are specific to minimal area problem
minimal_surface_element.feRhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpGetDataAtGaussPts("U",mse_common_data)
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpCalculateCoefficientsAtGaussPts("U",mse_common_data,true)
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpAssembleResidaul("U",mse_common_data)
);

and left hand side operator

// Calculate inverse of the 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));
minimal_surface_element.feLhs.getOpPtrVector().push_back(new MoFEM::OpSetInvJacH1ForFace(inv_jac_ptr));
// Below operators are specific to minimal area problem
minimal_surface_element.feLhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpGetDataAtGaussPts("U",mse_common_data)
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpCalculateCoefficientsAtGaussPts("U",mse_common_data,false)
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
  new MinimalSurfaceEquation::MinimalSurfaceElement::OpAssembleTangent("U",mse_common_data)
);

Note that general purpose operators MoFEM::OpCalculateJacForFace, OpInvertMatrix, and MoFEM::OpSetInvJacH1ForFace are used to calculate inverse of the Jacobian, and apply it to gradient of shape functions. Those two operators are implemented as a part core MoFEM library and are used in different contexts if needed.

All operators are run in order which are pushed to given element, so that first inverse of the Jacobian is calculated, then transverse of the inverse of the the Jacobian matrix it is applied to gradients of shape functions, next field gradient is calculated, then coefficients and finally residual vector or tangent matrix.

The elements for the right hand side and the left hand side calculations are added to SNES elements using following code;

DMMoFEMSNESSetFunction(dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feRhs,PETSC_NULL,PETSC_NULL);
DMMoFEMSNESSetJacobian(dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feLhs,PETSC_NULL,PETSC_NULL);

Example

Running code:

cd $MOFEM_INSTALL_DIR/users_modules/minimal_surface_equationa
mpirun -np 4 ./minimal_surface_area \
-my_file meshes/circle.h5m \
-ksp_monitor -ksp_converged_reason \
-pc_type lu -pc_factor_mat_solver_package mumps -mat_mumps_icntl_20 0 \
-snes_linesearch_type l2 -snes_linesearch_monitor -snes_linesearch_minlambda 0.01 \
-snes_monitor -snes_converged_reason \
-my_order 4 -my_nb_sub_steps 60 && \
mbconvert out.h5m out.vtk

Note

If you woring in build directory meshes are located in moudle source directory. If you working in install directory meshes are located in your current directory.

Program minimal_surface_equationa is started using mpirun, where the number of processors is selected to four in above example. The key parameters are order or approximation and number of sub-steps. In addition user can select different types of line-searchers or change SNES solver and select trust method. More details are available in PETSc documentation of SNES solver, see for details http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESSetFromOptions.html.

In this particular example is run for circular wire loop, function on boundary is given by

\[ u = \sin \left[ \pi (x+y) \right],\quad \textrm{on}\,\partial\Omega \]

Soap film spanned between circular wire

Mesh file is created with Cubit https://cubit.sandia.gov

reset
#create surface circle radius 1 zplane 
create surface rectangle width 2 height 2 zplane
surface all scheme TriMesh
surface all size auto factor 6
mesh surface all
refine curve all  numsplit 2 bias 1.0 depth 1 smooth

Future work and contribution

If you like to contribute to this module (or MoFEM library), you are very welcome. Pleas feel free to add changes, including improving this documentation. If you like to work and modify this user module, you can fork this repository, make your changes, and do pull request. Contact us cmatg.nosp@m.u@go.nosp@m.ogleg.nosp@m.roup.nosp@m.s.com or https://mofem.slack.com/ if you need any help or comments.

This is work in progress and more stuff can be done. You could consider own paper or student project.

Todo:

Make this documentation better

Error estimator for L2, H1 and H1 semi-norm

Mesh hp-adaptivity (All tools are there, need be used in this example).

Generalization to arbitrary wire curve in 3d space.

Write documentation about mesh bit levels and mesh refining.

Convergence analysis h- vs. p- refinement.

Adaptive time stepping

The plain program

/** \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;
}