v0.13.2
Loading...
Searching...
No Matches
SCL-5: Minimal surface equation
Note
Prerequisite of this tutorial is SCL-4: Nonlinear Poisson's equation


Note
Intended learning outcome:
  • practice of solving a nonlinear problem in MoFEM

Introduction

Introduction

Note
A different tutorial written in an old way solving the same problem can be found at Soap film spanned on wire

Plain program

The plain program for both the implementation of the UDOs (*.cpp) and the main program (*.cpp) are as follows

Implementation of the main program (*.cpp)

#include <stdlib.h>
using namespace MoFEM;
static char help[] = "...\n\n";
/** \brief Integrate the domain tangent matrix (LHS)
\f[
\sum\limits_j {\left[ {\int\limits_{{\Omega _e}} {\left( {{a_n}\nabla {\phi _i}
\cdot \nabla {\phi _j} - a_n^3\nabla {\phi _i}\left( {\nabla u \cdot \nabla
{\phi _j}} \right)\nabla u} \right)d{\Omega _e}} } \right]\delta {U_j}} =
\int\limits_{{\Omega _e}} {{\phi _i}fd{\Omega _e}} - \int\limits_{{\Omega _e}}
{\nabla {\phi _i}{a_n}\nabla ud{\Omega _e}} \\
{a_n} = \frac{1}{{{{\left( {1 +
{{\left| {\nabla u} \right|}^2}} \right)}^{\frac{1}{2}}}}}
\f]
*/
public:
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name,
boost::shared_ptr<MatrixDouble> field_grad_mat)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
fieldGradMat(field_grad_mat) {}
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get gradient of the field at integration points
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
// get derivatives of base functions on row
auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(i) * t_field_grad(i));
for (int rr = 0; rr != AssemblyDomainEleOp::nbRows; ++rr) {
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
for (int cc = 0; cc != AssemblyDomainEleOp::nbRows; cc++) {
// calculate components of the local matrix
locLhs(rr, cc) += (t_row_diff_base(i) * t_col_diff_base(i)) * an * a -
t_row_diff_base(i) *
(t_field_grad(i) * t_col_diff_base(i)) *
t_field_grad(i) * an * an * an * a;
// move to the derivatives of the next base functions on column
++t_col_diff_base;
}
// move to the derivatives of the next base functions on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
// move to the gradient of the field at the next integration point
++t_field_grad;
}
}
private:
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
/** \brief Integrate the domain residual vector (RHS)
\f[
\sum\limits_j {\left[ {\int\limits_{{\Omega _e}} {\left( {{a_n}\nabla {\phi _i}
\cdot \nabla {\phi _j} - a_n^3\nabla {\phi _i}\left( {\nabla u \cdot \nabla
{\phi _j}} \right)\nabla u} \right)d{\Omega _e}} } \right]\delta {U_j}} =
\int\limits_{{\Omega _e}} {{\phi _i}fd{\Omega _e}} - \int\limits_{{\Omega _e}}
{\nabla {\phi _i}{a_n}\nabla ud{\Omega _e}} \\
{a_n} = \frac{1}{{{{\left( {1 +
{{\left| {\nabla u} \right|}^2}} \right)}^{\frac{1}{2}}}}}
\f]
*/
public:
boost::shared_ptr<MatrixDouble> field_grad_mat)
fieldGradMat(field_grad_mat) {}
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get gradient of the field at integration points
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
// get base functions
auto t_base = data.getFTensor0N();
// get derivatives of base functions
auto t_diff_base = data.getFTensor1DiffN<2>();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(i) * t_field_grad(i));
for (int rr = 0; rr != AssemblyDomainEleOp::nbRows; rr++) {
// calculate components of the local vector
// remember to use -= here due to PETSc consideration of Residual Vec
nf[rr] += (t_diff_base(i) * t_field_grad(i)) * an * a;
// move to the next base function
++t_base;
// move to the derivatives of the next base function
++t_diff_base;
}
// move to the weight of the next integration point
++t_w;
// move to the gradient of the field at the next integration point
++t_field_grad;
}
}
private:
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(2 * M_PI * (x + y));
}
// Main interfaces
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
: mField(m_field) {}
}
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
CHKERR simple->addBoundaryField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setUp();
}
auto integration_rule = [](int o_row, int o_col, int approx_order) {
return 2 * approx_order;
};
auto *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
}
auto get_ents_on_mesh_skin = [&]() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, 2, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
// filter not owned entities, those are not on boundary
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
Range skin_verts;
mField.get_moab().get_connectivity(boundary_ents, skin_verts, true);
boundary_ents.merge(skin_verts);
boundaryEnts = boundary_ents;
return boundary_ents;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto problem_manager = mField.getInterface<ProblemsManager>();
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
problem_manager->markDofs(simple->getProblemName(), ROW,
ProblemsManager::OR, skin_edges, *marker_ptr);
return marker_ptr;
};
// Get global local vector of marked DOFs. Is global, since is set for all
// DOFs on processor. Is local since only DOFs on processor are in the
// vector. To access DOFs use local indices.
boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
}
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
pipeline.push_back(
new OpDomainTangentMatrix("U", "U", grad_u_at_gauss_pts));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
pipeline.push_back(new OpDomainResidualVector("U", grad_u_at_gauss_pts));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_base_ops = [&](auto &pipeline) {};
auto add_lhs_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
pipeline.push_back(new OpBoundaryMass(
"U", "U", [](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_rhs_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
pipeline.push_back(new OpBoundaryTimeScalarField(
"U", u_at_gauss_pts,
[](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpBoundarySource("U", boundaryFunction));
pipeline.push_back(new OpUnSetBc("U"));
};
auto pipeline_mng = mField.getInterface<PipelineManager>();
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto set_fieldsplit_preconditioner = [&](auto snes) {
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
SmartPetscObj<IS> is_all_bc;
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
name_prb, ROW, "U", 0, 1, is_all_bc, &boundaryEnts);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
}
};
// Create global RHS and solution vectors
auto dm = simple->getDM();
SmartPetscObj<Vec> global_rhs, global_solution;
global_solution = smartVectorDuplicate(global_rhs);
// Create nonlinear solver (SNES)
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createSNES();
CHKERR SNESSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR SNESSetUp(solver);
// Solve the system
CHKERR SNESSolve(solver, global_rhs, global_solution);
// Scatter result data on the mesh
CHKERR DMoFEMMeshToGlobalVector(dm, global_solution, INSERT_VALUES,
SCATTER_REVERSE);
}
auto post_proc = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", u_ptr));
post_proc->getOpPtrVector().push_back(
new OpPPMap(post_proc->getPostProcMesh(),
post_proc->getMapGaussPts(),
{{"U", u_ptr}},
{}, {}, {}
)
);
auto *simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), post_proc);
CHKERR post_proc->writeFile("out_result.h5m");
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example")
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
MinimalSurfaceEqn minimal_surface_problem(m_field);
CHKERR minimal_surface_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:574
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1155
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:523
SmartPetscObj< SNES > createSNES(SmartPetscObj< DM > dm=nullptr)
Create SNES (nonlinear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:33
FTensor::Index< 'i', 2 > i
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()
static double boundaryFunction(const double x, const double y, const double z)
MinimalSurfaceEqn(MoFEM::Interface &m_field)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Integrate the domain residual vector (RHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Integrate the domain tangent matrix (LHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.