v0.12.1
Soap film spanned on wire

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
constexpr double R

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();
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:71
@ H1
continuous field
Definition: definitions.h:96
#define CHKERR
Inline error check.
Definition: definitions.h:528

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"
);
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85

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 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);
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:458
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59

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 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);
constexpr double A
const double T
const FTensor::Tensor2< T, Dim, Dim > Vec

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
);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
CHKERR MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:540

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);
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:490

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);
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:670
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:711

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);
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:986
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:126
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:17

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);
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:478

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

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

OpAssmebleBcRhs(const string field_name,Vec v_f):
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
) {
int nb_dofs = data.getFieldData().size();
if(nb_dofs == 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);
}
vF,data.getIndices().size(),&data.getIndices()[0],&nF[0],ADD_VALUES
);
}
};
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:440
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:339
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:409
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual double evalFunction(double x, double y)
Function on boundary Inherit this class and overload this function to change bc.

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):
mA(m_a) {}
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type,EntityType col_type,
) {
int row_nb_dofs = row_data.getIndices().size();
if(row_nb_dofs == 0) {
}
int col_nb_dofs = col_data.getIndices().size();
if(col_nb_dofs == 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);
}
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);
mA,
col_nb_dofs,&*col_data.getIndices().data().begin(),
row_nb_dofs,&*row_data.getIndices().data().begin(),
&*transK.data().begin()
,ADD_VALUES
);
}
}
};
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:109
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)

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();
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73

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;

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 gradU
nb_gauss_pts x 2 gradients at integration pts
VectorDouble normGradU2
size of nb_gauss_pts, norm of gradient

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

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;

OpGetDataAtGaussPts(const string field_name,CommonData &common_data):
commonData(common_data) {
}
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
int nb_dofs = data.getFieldData().size();
if(nb_dofs == 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());
}
}
};
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:126
Evaluate function values and gradients at Gauss Pts.
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetDataAtGaussPts(const string field_name, CommonData &common_data)

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;

OpCalculateCoefficientsAtGaussPts(const string field_name,CommonData &common_data,bool rhs_and_not_lhs):
commonData(common_data),
rhsAndNotLhs(rhs_and_not_lhs) {
}
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
// 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]);
}
// 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];
}
}
}
}
}
};
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpCalculateCoefficientsAtGaussPts(const string field_name, CommonData &common_data, bool rhs_and_not_lhs)

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;

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

and tangent matrix

OpAssembleTangent(const string field_name,CommonData &common_data):
commonData(common_data) {
sYmm = false; // matrix is not symmetric
}
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type, EntityType col_type,
) {
int row_nb_dofs = row_data.getIndices().size();
if(row_nb_dofs == 0) {
}
int col_nb_dofs = col_data.getIndices().size();
if(col_nb_dofs == 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));
}
getFEMethod()->snes_B,
row_nb_dofs,&*row_data.getIndices().data().begin(),
col_nb_dofs,&*col_data.getIndices().data().begin(),
&*kTangent.data().begin()
,ADD_VALUES
);
}
};
OpAssembleTangent(const string field_name, CommonData &common_data)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)

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
// 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(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
minimal_surface_element.feRhs.getOpPtrVector().push_back(
);
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpPlasticTools::CommonData CommonData

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(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);
minimal_surface_element.feLhs.getOpPtrVector().push_back(
);

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 \
-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 <DirichletBC.hpp>
using namespace MinimalSurfaceEquation;
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
char mesh_file_name[255];
PetscInt order = 2;
PetscInt nb_sub_steps = 10;
PetscBool flg_file = PETSC_FALSE;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_atom_test = PETSC_FALSE;
{
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "",
"Minimal surface area config", "none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order", "default approximation order", "", 2,
&order, PETSC_NULL);
CHKERR PetscOptionsInt("-my_nb_sub_steps", "number of substeps", "", 10,
&nb_sub_steps, PETSC_NULL);
CHKERR PetscOptionsBool("-my_is_partitioned",
"set if mesh is partitioned (this result that "
"each process keeps only part of the mesh",
"", PETSC_FALSE, &is_partitioned, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_atom_test",
"is used with testing, exit with error when diverged", "",
PETSC_FALSE, &is_atom_test, PETSC_NULL);
ierr = PetscOptionsEnd();
// Reade parameters from line command
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
if (is_partitioned == PETSC_TRUE) {
// Read mesh to MOAB
const char *option;
option = "PARALLEL=BCAST_DELETE;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR moab.load_file(mesh_file_name, 0, option);
} else {
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
}
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
// meshset consisting all entities in mesh
EntityHandle root_set = moab.get_root_set();
// Seed all mesh entities to MoFEM database, those entities can be
// potentially used as finite elements or as entities which carry some
// approximation field.
BitRefLevel bit_level0;
bit_level0.set(0);
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 2, bit_level0);
// Add field
// Add entities to field (root_mesh, i.e. on all mesh entities fields are
// approx.) On those entities and lower dimension entities approximation is
// spaned,
CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "U");
// Set approximation oder
CHKERR m_field.set_field_order(root_set, MBTRI, "U", order);
CHKERR m_field.set_field_order(root_set, MBEDGE, "U", order);
CHKERR m_field.set_field_order(root_set, MBVERTEX, "U", 1);
CHKERR m_field.build_fields();
// Add finite element (this defines element declaration comes later)
CHKERR m_field.add_finite_element("MINIMAL_SURFACE_ELEMENT");
// Set on which fields that element operates
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
"MINIMAL_SURFACE_ELEMENT", "U");
// Add entities to that element
root_set, MBTRI, "MINIMAL_SURFACE_ELEMENT");
CHKERR m_field.add_finite_element("BC_ELEMENT");
CHKERR m_field.modify_finite_element_add_field_row("BC_ELEMENT", "U");
CHKERR m_field.modify_finite_element_add_field_col("BC_ELEMENT", "U");
CHKERR m_field.modify_finite_element_add_field_data("BC_ELEMENT", "U");
// Add entities to that element
Range bc_edges, bc_ents;
{
// Take a skin form the mesh
Skinner skin(&m_field.get_moab());
Range tris;
CHKERR moab.get_entities_by_type(root_set, MBTRI, tris, false);
Range bc_edges;
CHKERR skin.find_skin(root_set, tris, false, bc_edges);
Range bc_nodes;
CHKERR moab.get_connectivity(bc_edges, bc_nodes);
bc_ents.merge(bc_edges);
bc_ents.merge(bc_nodes);
// Add entities to that element
CHKERR m_field.add_ents_to_finite_element_by_type(bc_edges, MBEDGE,
"BC_ELEMENT");
}
// build finite elements
// build adjacencies between elements and degrees of freedom
CHKERR m_field.build_adjacencies(bit_level0);
// Create minimal surface area and elements
MinimalSurfaceElement minimal_surface_element(m_field);
// Data structure for common data
// This class is used to fix dofs on entities
DirichletFixFieldAtEntitiesBc fix_edges_ents(m_field, "U", NULL, NULL, NULL,
bc_ents);
// 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);
&fix_edges_ents);
// Lhs
CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, NULL, &fix_edges_ents,
NULL);
CHKERR DMMoFEMSNESSetJacobian(dm, "MINIMAL_SURFACE_ELEMENT",
&minimal_surface_element.feLhs, PETSC_NULL,
PETSC_NULL);
&fix_edges_ents);
}
SNESConvergedReason snes_reason;
// Solve problem
{
SNES snes;
SnesCtx *snes_ctx;
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
// CHKERR SNESSetDM(snes,dm);
CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
Vec T0;
CHKERR VecDuplicate(T, &T0);
CHKERR DMoFEMMeshToLocalVector(dm, T0, INSERT_VALUES, SCATTER_FORWARD);
double step_size = 1. / nb_sub_steps;
for (int ss = 1; ss <= nb_sub_steps; ss++) {
CHKERR VecAXPY(T, step_size, T0);
CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
CHKERR SNESSolve(snes, PETSC_NULL, T);
CHKERR SNESGetConvergedReason(snes, &snes_reason);
int its;
CHKERR SNESGetIterationNumber(snes, &its);
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"number of Newton iterations = %D\n\n", its);
if (snes_reason < 0) {
if (is_atom_test) {
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged");
} else {
break;
}
}
}
CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
// CHKERR MatView(A,PETSC_VIEWER_DRAW_SELF);
// std::string wait;
// std::cin >> wait;
CHKERR VecDestroy(&T0);
CHKERR SNESDestroy(&snes);
}
{
PostProcFaceOnRefinedMesh post_proc(m_field);
CHKERR 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);
}
}
return 0;
}
Implementation of minimal surface element.
Post-process fields on refined mesh.
int main(int argc, char *argv[])
static char help[]
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:365
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:208
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:476
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:51
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
char mesh_file_name[255]
Minimal surface equation.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1098
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
Fix dofs on entities.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:93
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.