|
| v0.14.0
|
- Note
- Prerequisite of this tutorial is SCL-3: Poisson's equation (Lagrange multiplier)
- Note
- Intended learning outcome:
- solving a nonlinear problem in MoFEM
- usage of PETSc SNES solver to solve nonlinear equations
Introduction
Plain program
The plain program for both the implementation of the UDOs (*
.hpp) and the main program (*
.cpp) are as follows
Implementation of User Data Operators (*.hpp)
#ifndef __NONLINEARPOISSON2D_HPP__
#define __NONLINEARPOISSON2D_HPP__
#include <stdlib.h>
typedef boost::function<
double(
const double,
const double,
const double)>
struct DataAtGaussPoints {
};
public:
std::string row_field_name, std::string col_field_name,
boost::shared_ptr<DataAtGaussPoints> &common_data,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
}
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locLhs.resize(nb_row_dofs, nb_col_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
auto t_field_grad = getFTensor1FromMat<2>(
commonData->fieldGrad);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += (((1 + t_field * t_field) * t_row_diff_base(
i) *
(2.0 * t_field * t_field_grad(
i) *
t_row_diff_base(
i) * t_col_base)) *
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
++t_field;
++t_field_grad;
}
}
}
}
ADD_VALUES);
}
}
private:
};
public:
boost::shared_ptr<DataAtGaussPoints> &common_data,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
if (nb_dofs) {
locRhs.resize(nb_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
auto t_field_grad = getFTensor1FromMat<2>(
commonData->fieldGrad);
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
t_grad_base(
i) * t_field_grad(
i) * (1 + t_field * t_field)) *
++t_base;
++t_grad_base;
}
++t_w;
++t_coords;
++t_field;
++t_field_grad;
}
}
}
}
CHKERR VecSetOption(
getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
}
private:
};
struct OpBoundaryTangentMatrix :
public OpEdgeEle {
public:
std::string col_field_name)
}
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
++t_col_base;
}
++t_row_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
}
}
}
private:
};
struct OpBoundaryResidualVector :
public OpEdgeEle {
public:
boost::shared_ptr<DataAtGaussPoints> &common_data)
if (nb_dofs) {
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
double boundary_term =
for (int rr = 0; rr != nb_dofs; rr++) {
++t_base;
}
++t_w;
++t_coords;
++t_field;
}
ADD_VALUES);
}
}
private:
};
};
#endif //__NONLINEARPOISSON2D_HPP__
Implementation of the main program (*.cpp)
#include <stdlib.h>
static char help[] =
"...\n\n";
public:
private:
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
}
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
}
MPI_Comm mpiComm;
const int mpiRank;
std::string domainField;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<FaceEle> domainTangentMatrixPipeline;
boost::shared_ptr<FaceEle> domainResidualVectorPipeline;
boost::shared_ptr<EdgeEle> boundaryTangentMatrixPipeline;
boost::shared_ptr<EdgeEle> boundaryResidualVectorPipeline;
boost::shared_ptr<DataAtGaussPoints> previousUpdate;
boost::shared_ptr<VectorDouble> fieldValuePtr;
boost::shared_ptr<MatrixDouble> fieldGradPtr;
boost::shared_ptr<PostProcEle> postProc;
Range boundaryEntitiesForFieldsplit;
};
: domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
mpiRank(mField.get_comm_rank()) {
domainTangentMatrixPipeline = boost::shared_ptr<FaceEle>(
new FaceEle(mField));
domainResidualVectorPipeline =
boost::shared_ptr<FaceEle>(
new FaceEle(mField));
boundaryTangentMatrixPipeline =
boost::shared_ptr<EdgeEle>(
new EdgeEle(mField));
boundaryResidualVectorPipeline =
boost::shared_ptr<EdgeEle>(
new EdgeEle(mField));
previousUpdate =
fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldValue);
fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
&previousUpdate->fieldGrad);
}
}
}
}
auto domain_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * (p + 1); };
auto domain_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * (p + 1); };
auto boundary_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * p + 1; };
auto boundary_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * p + 1; };
}
auto get_ents_on_mesh_skin = [&]() {
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
boundary_entities, true);
}
}
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](
Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{
}
{
}
{
}
{
}
{
boost::shared_ptr<FaceEle> null_face;
null_face);
null_face);
boost::shared_ptr<EdgeEle> null_edge;
null_edge);
null_edge);
}
}
if (1) {
KSP ksp_solver;
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
IS is_boundary;
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
SCATTER_REVERSE);
}
auto u_ptr = boost::make_shared<VectorDouble>();
{{"U", u_ptr}},
{}, {}, {}
)
);
dM, simpleInterface->getDomainFEName(),
boost::dynamic_pointer_cast<FEMethod>(postProc));
CHKERR postProc->writeFile(
"out_result.h5m");
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
try {
DMType dm_name = "DMMOFEM";
CHKERR poisson_problem.runProgram();
}
return 0;
}
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode outputResults()
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
NonlinearPoisson(MoFEM::Interface &m_field)
Integrate the domain tangent matrix (LHS)
OpDomainResidualVector(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
Problem manager is used to build and partition problems.
boost::shared_ptr< VectorDouble > fieldValuePtr
static double boundaryFunction(const double x, const double y, const double z)
virtual MPI_Comm & get_comm() const =0
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
auto createSNES(MPI_Comm comm)
MoFEMErrorCode boundaryCondition()
static double sourceTermFunction(const double x, const double y, const double z)
Integrate the domain residual vector (RHS)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
VectorDouble locBoundaryRhs
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Simple interface for fast problem set-up.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
@ OPROWCOL
operator doWork is executed on FE rows &columns
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< DataAtGaussPoints > previousUpdate
boost::shared_ptr< EdgeEle > boundaryTangentMatrixPipeline
auto getFTensor0IntegrationWeight()
Get integration weights.
Deprecated interface functions.
double getMeasure() const
get measure of element
ScalarFunc sourceTermFunc
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
DeprecatedCoreInterface Interface
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode getDM(DM *dm)
Get DM.
#define CHKERR
Inline error check.
friend class UserDataOperator
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
virtual moab::Interface & get_moab()=0
OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function, boost::shared_ptr< DataAtGaussPoints > &common_data)
MoFEMErrorCode assembleSystem()
implementation of Data Operators for Forces and Sources
Section manager is used to create indexes and sections.
boost::shared_ptr< FaceEle > domainResidualVectorPipeline
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
SmartPetscObj< SNES > snesSolver
boost::shared_ptr< DataAtGaussPoints > commonData
default operator for TRI element
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode solveSystem()
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
FTensor::Index< 'i', 2 > i
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Get value at integration points for scalar field.
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MatrixDouble transLocBoundaryLhs
const std::string getDomainFEName() const
Get the Domain FE Name.
Modify integration weights on face to take in account higher-order geometry.
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.
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
boost::shared_ptr< MatrixDouble > fieldGradPtr
friend class UserDataOperator
boost::shared_ptr< DataAtGaussPoints > commonData
int main(int argc, char *argv[])
MoFEMErrorCode runProgram()
MoFEM::Interface & mField
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
EntitiesFieldData::EntData EntData
constexpr auto field_name
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
default operator for EDGE element
MatrixDouble locBoundaryLhs
MoFEM::FaceElementForcesAndSourcesCore FaceEle
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
#define CATCH_ERRORS
Catch errors.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
OpBoundaryTangentMatrix(std::string row_field_name, std::string col_field_name)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
boost::shared_ptr< FaceEle > domainTangentMatrixPipeline
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name, boost::shared_ptr< DataAtGaussPoints > &common_data, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
UBlasVector< double > VectorDouble
MoFEMErrorCode addBoundaryField(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_ZERO, int verb=-1)
Add field on boundary.
MoFEMErrorCode setupProblem()
boost::shared_ptr< DataAtGaussPoints > commonData
keeps basic data about problem
boost::shared_ptr< EdgeEle > boundaryResidualVectorPipeline
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
boost::function< double(const double, const double, const double)> ScalarFunc
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
bool sYmm
If true assume that matrix is symmetric structure.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode readMesh()
Range boundaryEntitiesForFieldsplit
@ OPROW
operator doWork function is executed on FE rows
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Post post-proc data at points from hash maps.