v0.6.9
Soap film spanned on wire

Table of Contents

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_rectangle.png
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
ierr = m_field.add_field("U",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
// 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,
ierr = m_field.add_ents_to_field_by_type(root_set,MBTRI,"U"); CHKERRQ(ierr);
// Set approximation oder
ierr = m_field.set_field_order(root_set,MBTRI,"U",order); CHKERRQ(ierr);
ierr = m_field.set_field_order(root_set,MBEDGE,"U",order); CHKERRQ(ierr);
ierr = m_field.set_field_order(root_set,MBVERTEX,"U",1); CHKERRQ(ierr);
ierr = m_field.build_fields(); CHKERRQ(ierr);

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)
ierr = m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);
// Set (define) on which fields that element operates
ierr = m_field.modify_finite_element_add_field_row("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);

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

ierr = m_field.add_finite_element("BC_ELEMENT"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("BC_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("BC_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("BC_ELEMENT","U"); CHKERRQ(ierr);

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
ierr = m_field.add_ents_to_finite_element_by_type(root_set,MBTRI,"MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);

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); CHKERRQ_MOAB(rval);
Range bc_edges;
rval = skin.find_skin(root_set,tris,false,bc_edges); CHKERRQ_MOAB(rval);
// Add entities to that element
bc_edges,"BC_ELEMENT"
); CHKERRQ(ierr);

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
ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD,&bc_dm);CHKERRQ(ierr);
ierr = DMSetType(bc_dm,dm_name);CHKERRQ(ierr);
// Create DM instance
ierr = DMMoFEMCreateMoFEM(bc_dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
// Configure DM from line command options (DM itself, solvers, pre-conditioners, ... )
ierr = DMSetFromOptions(bc_dm); CHKERRQ(ierr);
// Add elements to dm (only one here)
ierr = DMMoFEMAddElement(bc_dm,"BC_ELEMENT"); CHKERRQ(ierr);
// Set up problem (number dofs, partition mesh, etc.)
ierr = DMSetUp(bc_dm); CHKERRQ(ierr);

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
ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
// Create DM instance
ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
// Get options fron line command
ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
// Add elements to dm (only one here)
ierr = DMMoFEMAddElement(dm,"MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);
// Set up data structures
ierr = DMSetUp(dm); CHKERRQ(ierr);

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
ierr = DMCreateGlobalVector(dm,&T); CHKERRQ(ierr);
ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
// Create matrix
ierr = DMCreateMatrix(dm,&A); CHKERRQ(ierr);

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
bc_dm,"BC_ELEMENT",&minimal_surface_element.feBcEdge
); CHKERRQ(ierr);
ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

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;
ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
ierr = KSPSetUp(solver); CHKERRQ(ierr);
ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
ierr = KSPDestroy(&solver); CHKERRQ(ierr);

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

ierr = DMoFEMMeshToGlobalVector(bc_dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);

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
ierr = DMMoFEMSNESSetFunction(dm,DM_NO_ELEMENT,NULL,&fix_edges_ents,NULL); CHKERRQ(ierr);
ierr = DMMoFEMSNESSetFunction(dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feRhs,PETSC_NULL,PETSC_NULL);
ierr = DMMoFEMSNESSetFunction(dm,DM_NO_ELEMENT,NULL,NULL,&fix_edges_ents); CHKERRQ(ierr);
//Lhs
ierr = DMMoFEMSNESSetJacobian(dm,DM_NO_ELEMENT,NULL,&fix_edges_ents,NULL); CHKERRQ(ierr);
ierr = DMMoFEMSNESSetJacobian(dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feLhs,PETSC_NULL,PETSC_NULL);
ierr = DMMoFEMSNESSetJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&fix_edges_ents); CHKERRQ(ierr);

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;
ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
//ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,A,A,SnesMat,snes_ctx); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

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;
ierr = VecDuplicate(T,&T0); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm,T0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
double step_size = 1./nb_sub_steps;
for(int ss = 1;ss<=nb_sub_steps;ss++) {
ierr = VecAXPY(T,step_size,T0); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = SNESSolve(snes,PETSC_NULL,T); CHKERRQ(ierr);
}
ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);

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 [41]. 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

/** \brief Set integrate rule
*/
int getRule(int order) { return 2*order; };
};

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

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));
}
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
PetscFunctionBegin;
int nb_dofs = data.getFieldData().size();
if(nb_dofs == 0) {
PetscFunctionReturn(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);
}
ierr = VecSetValues(
vF,data.getIndices().size(),&data.getIndices()[0],&nF[0],ADD_VALUES
); CHKERRQ(ierr);
PetscFunctionReturn(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,

Mat mA;
OpAssmebleBcLhs(const string field_name,Mat m_a):
EdgeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROWCOL),
mA(m_a) {}
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
int row_nb_dofs = row_data.getIndices().size();
if(row_nb_dofs == 0) {
PetscFunctionReturn(0);
}
int col_nb_dofs = col_data.getIndices().size();
if(col_nb_dofs == 0) {
PetscFunctionReturn(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);
}
ierr = MatSetValues(
mA,
row_nb_dofs,&*row_data.getIndices().data().begin(),
col_nb_dofs,&*col_data.getIndices().data().begin(),
&*K.data().begin()
,ADD_VALUES
); CHKERRQ(ierr);
// 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);
ierr = MatSetValues(
mA,
col_nb_dofs,&*col_data.getIndices().data().begin(),
row_nb_dofs,&*row_data.getIndices().data().begin(),
&*transK.data().begin()
,ADD_VALUES
); CHKERRQ(ierr);
}
PetscFunctionReturn(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;

ierr = DMoFEMLoopFiniteElements(bc_dm,"BC_ELEMENT",&minimal_surface_element.feBcEdge); CHKERRQ(ierr);

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;

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$
MatrixDouble invJac; ///< inverse of the Jacobian matrix
};

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

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

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
) {
PetscFunctionBegin;
int nb_dofs = data.getFieldData().size();
if(nb_dofs == 0) {
PetscFunctionReturn(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());
}
PetscFunctionReturn(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;

CommonData &commonData;
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
) {
PetscFunctionBegin;
// 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];
}
}
}
}
PetscFunctionReturn(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;

CommonData &commonData;
OpAssembleResidaul(const string field_name,CommonData &common_data):
FaceElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
commonData(common_data) {
}
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
PetscFunctionBegin;
int nb_dofs = data.getIndices().size();
if(nb_dofs == 0) {
PetscFunctionReturn(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);
}
ierr = VecSetValues(
getFEMethod()->snes_f,
data.getIndices().size(),
&*data.getIndices().data().begin(),
&rEsidual[0],
ADD_VALUES
); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
};

and tangent matrix

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
}
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type, EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data
) {
PetscFunctionBegin;
try {
int row_nb_dofs = row_data.getIndices().size();
if(row_nb_dofs == 0) {
PetscFunctionReturn(0);
}
int col_nb_dofs = col_data.getIndices().size();
if(col_nb_dofs == 0) {
PetscFunctionReturn(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));
}
ierr = 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
); CHKERRQ(ierr);
} catch (exception& ex) {
ostringstream ss;
ss << "thorw in method: " << ex.what() << " at line " << __LINE__ << " in file " << __FILE__;
SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
}
PetscFunctionReturn(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 MoFEM::OpCalculateInvJacForFace(mse_common_data.invJac));
// Apply invJac to shape functions
minimal_surface_element.feRhs.getOpPtrVector().push_back(new MoFEM::OpSetInvJacH1ForFace(mse_common_data.invJac));
// Below operators are specific to minimal area problem
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);

and left hand side operator

// Calculate inverse of the Jacobian
minimal_surface_element.feLhs.getOpPtrVector().push_back(new MoFEM::OpCalculateInvJacForFace(mse_common_data.invJac));
// Apply invJac to shape functions
minimal_surface_element.feLhs.getOpPtrVector().push_back(new MoFEM::OpSetInvJacH1ForFace(mse_common_data.invJac));
// Below operators are specific to minimal area problem
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);

Note that general purpose operators MoFEM::OpCalculateInvJacForFace 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 ../users_modules/minimal_surface_equation/meshes/circle.h5m \
-ksp_monitor -ksp_converged_reason -pc_type lu -pc_factor_mat_solver_package mumps \
-snes_linesearch_type l2 -snes_linesearch_monitor -snes_linesearch_minlambda 0.01 \
-snes_monitor -snes_converged_reason \
-my_order 4 -my_nb_sub_steps 10 && \
mbconvert out.h5m out.vtk

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 \]

minimal_surface_equation.gif
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 gel_analysis.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 <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <moab/Skinner.hpp>
#include <DirichletBC.hpp>
using namespace MinimalSurfaceEquation;
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
PetscInitialize(&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;
{
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Minimal surface area config","none"); CHKERRQ(ierr);
ierr = PetscOptionsString(
"-my_file",
"mesh file name","",
"mesh.h5m",mesh_file_name,255,&flg_file
); CHKERRQ(ierr);
ierr = PetscOptionsInt(
"-my_order",
"default approximation order","",
2,&order,PETSC_NULL
); CHKERRQ(ierr);
ierr = PetscOptionsInt(
"-my_nb_sub_steps",
"number of substeps","",
10,&nb_sub_steps,PETSC_NULL
); CHKERRQ(ierr);
ierr = 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
); CHKERRQ(ierr);
ierr = PetscOptionsBool(
"-my_is_atom_test",
"is used with testing, exit with error when diverged","",
PETSC_FALSE,&is_partitioned,PETSC_NULL
); CHKERRQ(ierr);
ierr = PetscOptionsEnd(); CHKERRQ(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;";
rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
} else {
const char *option;
option = "";
rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
}
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);
ierr = m_field.query_interface<BitRefManager>()->setBitRefLevelByDim(0,2,bit_level0); CHKERRQ(ierr);
// Add field
ierr = m_field.add_field("U",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
// 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,
ierr = m_field.add_ents_to_field_by_type(root_set,MBTRI,"U"); CHKERRQ(ierr);
// Set approximation oder
ierr = m_field.set_field_order(root_set,MBTRI,"U",order); CHKERRQ(ierr);
ierr = m_field.set_field_order(root_set,MBEDGE,"U",order); CHKERRQ(ierr);
ierr = m_field.set_field_order(root_set,MBVERTEX,"U",1); CHKERRQ(ierr);
ierr = m_field.build_fields(); CHKERRQ(ierr);
// Add finite element (this defines element declaration comes later)
ierr = m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);
// Set on which fields that element operates
ierr = m_field.modify_finite_element_add_field_row("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("MINIMAL_SURFACE_ELEMENT","U"); CHKERRQ(ierr);
// Add entities to that element
ierr = m_field.add_ents_to_finite_element_by_type(root_set,MBTRI,"MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);
ierr = m_field.add_finite_element("BC_ELEMENT"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_row("BC_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_col("BC_ELEMENT","U"); CHKERRQ(ierr);
ierr = m_field.modify_finite_element_add_field_data("BC_ELEMENT","U"); CHKERRQ(ierr);
// Add entities to that element
Range bc_edges,bc_ents;
{
// 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); CHKERRQ_MOAB(rval);
Range bc_edges;
rval = skin.find_skin(root_set,tris,false,bc_edges); CHKERRQ_MOAB(rval);
Range bc_nodes;
rval = moab.get_connectivity(bc_edges,bc_nodes); CHKERRQ(rval);
bc_ents.merge(bc_edges);
bc_ents.merge(bc_nodes);
// Add entities to that element
bc_edges,MBEDGE,"BC_ELEMENT"
); CHKERRQ(ierr);
}
//build finite elements
ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
//build adjacencies between elements and degrees of freedom
ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
// 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
// ierr = m_field.add_problem("DM_BC"); CHKERRQ(ierr);
// //set refinement level for problem
// ierr = m_field.modify_problem_ref_level_add_bit("DM_BC",bit_level0); CHKERRQ(ierr);
DM bc_dm;
DMType dm_name = "DM_BC";
// Register DM problem
ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD,&bc_dm);CHKERRQ(ierr);
ierr = DMSetType(bc_dm,dm_name);CHKERRQ(ierr);
// Create DM instance
ierr = DMMoFEMCreateMoFEM(bc_dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
// Configure DM form line command options (DM itself, solvers, pre-conditioners, ... )
ierr = DMSetFromOptions(bc_dm); CHKERRQ(ierr);
// Add elements to dm (only one here)
ierr = DMMoFEMAddElement(bc_dm,"BC_ELEMENT"); CHKERRQ(ierr);
// Set up problem (number dofs, partition mesh, etc.)
ierr = DMSetUp(bc_dm); CHKERRQ(ierr);
Mat A;
Vec T,F;
ierr = DMCreateGlobalVector(bc_dm,&T); CHKERRQ(ierr);
ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
ierr = DMCreateMatrix(bc_dm,&A); CHKERRQ(ierr);
// Boundary problem
minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
);
minimal_surface_element.feBcEdge.getOpPtrVector().push_back(
);
bc_dm,"BC_ELEMENT",&minimal_surface_element.feBcEdge
); CHKERRQ(ierr);
ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
{
KSP solver;
ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
ierr = KSPSetUp(solver); CHKERRQ(ierr);
ierr = KSPSolve(solver,F,T); CHKERRQ(ierr);
ierr = KSPDestroy(&solver); CHKERRQ(ierr);
}
// Scatter data on mesh
ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = DMoFEMMeshToGlobalVector(bc_dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = VecDestroy(&T); CHKERRQ(ierr);
ierr = VecDestroy(&F); CHKERRQ(ierr);
ierr = MatDestroy(&A); CHKERRQ(ierr);
ierr = DMDestroy(&bc_dm); CHKERRQ(ierr);
}
// //define problems
// ierr = m_field.add_problem("DM_MINIMAL_AREA"); CHKERRQ(ierr);
// //set refinement level for problem
// ierr = m_field.modify_problem_ref_level_add_bit("DM_MINIMAL_AREA",bit_level0); CHKERRQ(ierr);
// Create discrete manager instance
DM dm;
DMType dm_name = "DM_MINIMAL_AREA";
{
// Register problem
ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
// Create DM instance
ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
// Get options fron line command
ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
// Add elements to dm (only one here)
ierr = DMMoFEMAddElement(dm,"MINIMAL_SURFACE_ELEMENT"); CHKERRQ(ierr);
// Set up data structures
ierr = DMSetUp(dm); CHKERRQ(ierr);
}
// Create matrices and vectors used for analysis
Vec T,F;
Mat A;
{
ierr = DMCreateGlobalVector(dm,&T); CHKERRQ(ierr);
ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm,&A); CHKERRQ(ierr);
}
//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(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
// 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(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);
//Rhs
ierr = DMMoFEMSNESSetFunction(dm,DM_NO_ELEMENT,NULL,&fix_edges_ents,NULL); CHKERRQ(ierr);
dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feRhs,PETSC_NULL,PETSC_NULL
); CHKERRQ(ierr);
ierr = DMMoFEMSNESSetFunction(dm,DM_NO_ELEMENT,NULL,NULL,&fix_edges_ents); CHKERRQ(ierr);
//Lhs
ierr = DMMoFEMSNESSetJacobian(dm,DM_NO_ELEMENT,NULL,&fix_edges_ents,NULL); CHKERRQ(ierr);
dm,"MINIMAL_SURFACE_ELEMENT",&minimal_surface_element.feLhs,PETSC_NULL,PETSC_NULL
); CHKERRQ(ierr);
ierr = DMMoFEMSNESSetJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&fix_edges_ents); CHKERRQ(ierr);
}
SNESConvergedReason snes_reason;
// Solve problem
{
SNES snes;
SnesCtx *snes_ctx;
ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
//ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,A,A,SnesMat,snes_ctx); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
Vec T0;
ierr = VecDuplicate(T,&T0); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm,T0,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
double step_size = 1./nb_sub_steps;
for(int ss = 1;ss<=nb_sub_steps;ss++) {
ierr = VecAXPY(T,step_size,T0); CHKERRQ(ierr);
ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
ierr = SNESSolve(snes,PETSC_NULL,T); CHKERRQ(ierr);
ierr = SNESGetConvergedReason(snes,&snes_reason); CHKERRQ(ierr);
int its;
ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %D\n\n",its); CHKERRQ(ierr);
if(snes_reason<0) {
if(is_atom_test) {
SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"atom test diverged");
} else {
break;
}
}
}
ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
// ierr = MatView(A,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
// std::string wait;
// std::cin >> wait;
ierr = VecDestroy(&T0); CHKERRQ(ierr);
ierr = SNESDestroy(&snes); CHKERRQ(ierr);
}
{
PostProcFaceOnRefinedMesh post_proc(m_field);
ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
ierr = post_proc.addFieldValuesPostProc("U"); CHKERRQ(ierr);
ierr = DMoFEMLoopFiniteElements(dm,"MINIMAL_SURFACE_ELEMENT",&post_proc); CHKERRQ(ierr);
//Save data on mesh
rval = post_proc.postProcMesh.write_file(
"out.h5m","MOAB","PARALLEL=WRITE_PART"
}
// Clean and destroy
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
ierr = VecDestroy(&T); CHKERRQ(ierr);
ierr = VecDestroy(&F); CHKERRQ(ierr);
ierr = MatDestroy(&A); CHKERRQ(ierr);
}
} catch (MoFEMException const &e) {
SETERRQ(PETSC_COMM_SELF,e.errorCode,e.errorMessage);
}
PetscFinalize();
return 0;
}