v0.14.0
Loading...
Searching...
No Matches
SCL-4: Nonlinear Poisson's equation
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>
using EntData = EntitiesFieldData::EntData;
FTensor::Index<'i', 2> i;
typedef boost::function<double(const double, const double, const double)>
struct DataAtGaussPoints {
// This struct is for data required for the update of Newton's iteration
VectorDouble fieldValue;
MatrixDouble fieldGrad;
};
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)
: OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
commonData(common_data), boundaryMarker(boundary_marker) {
sYmm = false;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_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);
locLhs.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get solution (field value) at integration points
auto t_field = getFTensor0FromVec(commonData->fieldValue);
// get gradient of field at integration points
auto t_field_grad = getFTensor1FromMat<2>(commonData->fieldGrad);
// get base functions on row
auto t_row_base = row_data.getFTensor0N();
// 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;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get base functions on column
auto t_col_base = col_data.getFTensor0N(gg, 0);
// get derivatives of base functions on column
auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += (((1 + t_field * t_field) * t_row_diff_base(i) *
t_col_diff_base(i)) +
(2.0 * t_field * t_field_grad(i) *
t_row_diff_base(i) * t_col_base)) *
a;
// move to the next base functions on column
++t_col_base;
// move to the derivatives of the next base function on column
++t_col_diff_base;
}
// move to the next base function on row
++t_row_base;
// move to the derivatives of the next base function on row
++t_row_diff_base;
}
// move to the weight of the next integration point
++t_w;
// move to the solution (field value) at the next integration point
++t_field;
// move to the gradient of field value at the next integration point
++t_field_grad;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
// store original row indices
auto row_indices = row_data.getIndices();
// mark the boundary DOFs (as -1) and fill only domain DOFs
for (int r = 0; r != row_data.getIndices().size(); ++r) {
if ((*boundaryMarker)[row_data.getLocalIndices()[r]]) {
row_data.getIndices()[r] = -1;
}
}
}
// fill value to local stiffness matrix ignoring boundary DOFs
CHKERR MatSetValues(getSNESB(), row_data, col_data, &locLhs(0, 0),
ADD_VALUES);
// revert back row indices to the original
row_data.getIndices().swap(row_indices);
}
}
private:
boost::shared_ptr<DataAtGaussPoints> commonData;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
MatrixDouble locLhs;
};
public:
std::string field_name, ScalarFunc source_term_function,
boost::shared_ptr<DataAtGaussPoints> &common_data,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
sourceTermFunc(source_term_function), commonData(common_data),
boundaryMarker(boundary_marker) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
locRhs.resize(nb_dofs, false);
locRhs.clear();
// get element area
const double area = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get coordinates of the integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// get solution (field value) at integration point
auto t_field = getFTensor0FromVec(commonData->fieldValue);
// get gradient of field value of integration point
auto t_field_grad = getFTensor1FromMat<2>(commonData->fieldGrad);
// get base function
auto t_base = data.getFTensor0N();
// get derivatives of base function
auto t_grad_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;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
// calculate the local vector
for (int rr = 0; rr != nb_dofs; rr++) {
locRhs[rr] -=
(t_base * body_source -
t_grad_base(i) * t_field_grad(i) * (1 + t_field * t_field)) *
a;
// move to the next base function
++t_base;
// move to the derivatives of the next base function
++t_grad_base;
}
// move to the weight of the next integration point
++t_w;
// move to the coordinates of the next integration point
++t_coords;
// move to the solution (field value) at the next integration point
++t_field;
// move to the gradient of field value at the next integration point
++t_field_grad;
}
// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
// store original row indices
auto row_indices = data.getIndices();
// mark the boundary DOFs (as -1) and fill only domain DOFs
for (int r = 0; r != data.getIndices().size(); ++r) {
if ((*boundaryMarker)[data.getLocalIndices()[r]]) {
data.getIndices()[r] = -1;
}
}
}
// fill value to local vector ignoring boundary DOFs
CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
CHKERR VecSetValues(getSNESf(), data, &*locRhs.begin(), ADD_VALUES);
// revert back the indices
data.getIndices().swap(row_indices);
}
}
private:
boost::shared_ptr<DataAtGaussPoints> commonData;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
VectorDouble locRhs;
};
struct OpBoundaryTangentMatrix : public OpEdgeEle {
public:
OpBoundaryTangentMatrix(std::string row_field_name,
std::string col_field_name)
: OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
sYmm = true;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
// cerr << nb_row_dofs;
if (nb_row_dofs && nb_col_dofs) {
locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
locBoundaryLhs.clear();
// get (boundary) element length
const double edge = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get base function on row
auto t_row_base = row_data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
// get base function on column
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col_dofs; cc++) {
locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
// move to the next base function on column
++t_col_base;
}
// move to the next base function on row
++t_row_base;
}
// move to the weight of the next integration point
++t_w;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
CHKERR MatSetValues(getSNESB(), row_data, col_data, &locBoundaryLhs(0, 0),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocBoundaryLhs.resize(nb_col_dofs, nb_row_dofs, false);
CHKERR MatSetValues(getSNESB(), col_data, row_data,
&transLocBoundaryLhs(0, 0), ADD_VALUES);
}
// cerr << locBoundaryLhs << endl;
// cerr << transLocBoundaryLhs << endl;
}
}
private:
};
struct OpBoundaryResidualVector : public OpEdgeEle {
public:
OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function,
boost::shared_ptr<DataAtGaussPoints> &common_data)
boundaryFunc(boundary_function), commonData(common_data) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
locBoundaryRhs.resize(nb_dofs, false);
locBoundaryRhs.clear();
// get (boundary) element length
const double edge = getMeasure();
// get number of integration points
const int nb_integration_points = getGaussPts().size2();
// get integration weights
// get coordinates at integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// get solution (field value) at integration point
auto t_field = getFTensor0FromVec(commonData->fieldValue);
// get base function
auto t_base = data.getFTensor0N();
// START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
double boundary_term =
boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
// calculate the local vector
for (int rr = 0; rr != nb_dofs; rr++) {
locBoundaryRhs[rr] -= t_base * (boundary_term - t_field) * a;
// move to the next base function
++t_base;
}
// move to the weight of the next integration point
++t_w;
// move to the coordinates of the next integration point
++t_coords;
// move to the solution (field value) at the next integration point
++t_field;
}
// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
ADD_VALUES);
}
}
private:
boost::shared_ptr<DataAtGaussPoints> commonData;
VectorDouble locBoundaryRhs;
};
}; // namespace NonlinearPoissonOps
#endif //__NONLINEARPOISSON2D_HPP__
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
FTensor::Index< 'i', 2 > i
boost::function< double(const double, const double, const double)> ScalarFunc
const double body_source
int r
Definition sdf.py:8
constexpr auto field_name
bool sYmm
If true assume that matrix is symmetric structure.
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.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
OpBoundaryResidualVector(std::string field_name, ScalarFunc boundary_function, boost::shared_ptr< DataAtGaussPoints > &common_data)
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpBoundaryTangentMatrix(std::string row_field_name, std::string col_field_name)
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.
boost::shared_ptr< DataAtGaussPoints > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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)
boost::shared_ptr< DataAtGaussPoints > commonData
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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)
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.
Integrate the domain residual vector (RHS)
Integrate the domain tangent matrix (LHS)

Implementation of the main program (*.cpp)

#include <stdlib.h>
using namespace MoFEM;
using namespace NonlinearPoissonOps;
static char help[] = "...\n\n";
public:
// Declaration of the main function to run analysis
private:
// Declaration of other main functions called in runProgram()
// Function to calculate the Source term
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
// return 1;
}
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
// return 0;
}
// Main interfaces
// mpi parallel communicator
MPI_Comm mpiComm;
// Number of processors
const int mpiRank;
// Discrete Manager and nonliner SNES solver using SmartPetscObj
// Field name and approximation order
std::string domainField;
int order;
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
// MoFEM working Pipelines for LHS and RHS of domain and boundary
boost::shared_ptr<FaceEle> domainTangentMatrixPipeline;
boost::shared_ptr<FaceEle> domainResidualVectorPipeline;
boost::shared_ptr<EdgeEle> boundaryTangentMatrixPipeline;
boost::shared_ptr<EdgeEle> boundaryResidualVectorPipeline;
// Objects needed for solution updates in Newton's method
boost::shared_ptr<DataAtGaussPoints> previousUpdate;
boost::shared_ptr<VectorDouble> fieldValuePtr;
boost::shared_ptr<MatrixDouble> fieldGradPtr;
// Object needed for postprocessing
boost::shared_ptr<PostProcEle> postProc;
// Boundary entities marked for fieldsplit (block) solver - optional
};
: 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 =
boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldValue);
fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
&previousUpdate->fieldGrad);
}
}
}
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
}
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); };
domainTangentMatrixPipeline->getRuleHook = domain_rule_lhs;
domainResidualVectorPipeline->getRuleHook = domain_rule_rhs;
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; };
boundaryTangentMatrixPipeline->getRuleHook = boundary_rule_lhs;
boundaryResidualVectorPipeline->getRuleHook = boundary_rule_rhs;
}
// Get boundary edges marked in block named "BOUNDARY_CONDITION"
auto get_ents_on_mesh_skin = [&]() {
Range boundary_entities;
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
boundary_entities, true);
}
}
// Add vertices to boundary entities
Range boundary_vertices;
CHKERR mField.get_moab().get_connectivity(boundary_entities,
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
// Store entities for fieldsplit (block) solver
boundaryEntitiesForFieldsplit = boundary_entities;
return boundary_entities;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto problem_manager = mField.getInterface<ProblemsManager>();
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
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 det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{ // Push operators to the Pipeline that is responsible for calculating the
// domain tangent matrix (LHS)
// Add default operators to calculate inverse of Jacobian (needed for
// implementation of 2D problem but not 3D ones)
domainTangentMatrixPipeline->getOpPtrVector().push_back(
new OpCalculateHOJac<2>(jac_ptr));
domainTangentMatrixPipeline->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
domainTangentMatrixPipeline->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
domainTangentMatrixPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field values at integration points
domainTangentMatrixPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field gradient at integration points
domainTangentMatrixPipeline->getOpPtrVector().push_back(
// Push operators for domain tangent matrix (LHS)
domainTangentMatrixPipeline->getOpPtrVector().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating the
// domain residual vector (RHS)
// Add default operators to calculate inverse of Jacobian (needed for
// implementation of 2D problem but not 3D ones)
domainResidualVectorPipeline->getOpPtrVector().push_back(
new OpCalculateHOJac<2>(jac_ptr));
domainResidualVectorPipeline->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
domainResidualVectorPipeline->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
domainResidualVectorPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field values at integration points
domainResidualVectorPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field gradient at integration points
domainResidualVectorPipeline->getOpPtrVector().push_back(
domainResidualVectorPipeline->getOpPtrVector().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating the
// boundary tangent matrix (LHS)
boundaryTangentMatrixPipeline->getOpPtrVector().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating
// boundary residual vector (RHS)
// Add default operator to calculate field values at integration points
boundaryResidualVectorPipeline->getOpPtrVector().push_back(
boundaryResidualVectorPipeline->getOpPtrVector().push_back(
}
// get Discrete Manager (SmartPetscObj)
{ // Set operators for nonlinear equations solver (SNES) from MoFEM Pipelines
// Set operators for calculation of LHS and RHS of domain elements
boost::shared_ptr<FaceEle> null_face;
null_face);
null_face);
// Set operators for calculation of LHS and RHS of boundary elements
boost::shared_ptr<EdgeEle> null_edge;
null_edge);
null_edge);
}
}
// Create RHS and solution vectors
SmartPetscObj<Vec> global_rhs, global_solution;
global_solution = vectorDuplicate(global_rhs);
// Create nonlinear solver (SNES)
CHKERR SNESSetFromOptions(snesSolver);
// Fieldsplit block solver: yes/no
if (1) {
KSP ksp_solver;
CHKERR SNESGetKSP(snesSolver, &ksp_solver);
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Set up FIELDSPLIT, only when option used -pc_type fieldsplit
if (is_pcfs == PETSC_TRUE) {
IS is_boundary;
const MoFEM::Problem *problem_ptr;
CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
// CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
CHKERR SNESSetDM(snesSolver, dM);
CHKERR SNESSetUp(snesSolver);
// Solve the system
CHKERR SNESSolve(snesSolver, global_rhs, global_solution);
// VecView(global_rhs, PETSC_VIEWER_STDOUT_SELF);
// Scatter result data on the mesh
CHKERR DMoFEMMeshToGlobalVector(dM, global_solution, INSERT_VALUES,
SCATTER_REVERSE);
}
postProc = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
postProc->getOpPtrVector().push_back(
postProc->getOpPtrVector().push_back(
new OpPPMap(postProc->getPostProcMesh(), postProc->getMapGaussPts(),
{{"U", u_ptr}},
{}, {}, {}
)
);
dM, simpleInterface->getDomainFEName(),
boost::dynamic_pointer_cast<FEMethod>(postProc));
CHKERR postProc->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";
// 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
NonlinearPoisson poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
static Index< 'p', 3 > p
std::string param_file
static char help[]
int main()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
@ H1
continuous field
Definition definitions.h:85
@ BLOCKSET
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 DMMoFEM.cpp:704
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:412
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 DMMoFEM.cpp:745
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:572
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1153
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:521
#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.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto createSNES(MPI_Comm comm)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =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.
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
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 inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
keeps basic data about problem
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
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition Simple.hpp:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
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.
Definition Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
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.
Definition Simple.cpp:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:611
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:362
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:341
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode outputResults()
MoFEMErrorCode boundaryCondition()
static double boundaryFunction(const double x, const double y, const double z)
MoFEMErrorCode solveSystem()
SmartPetscObj< SNES > snesSolver
static double sourceTermFunction(const double x, const double y, const double z)
MoFEMErrorCode readMesh()
MoFEM::Interface & mField
SmartPetscObj< DM > dM
MoFEMErrorCode assembleSystem()
NonlinearPoisson(MoFEM::Interface &m_field)
boost::shared_ptr< EdgeEle > boundaryTangentMatrixPipeline
boost::shared_ptr< VectorDouble > fieldValuePtr
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
boost::shared_ptr< MatrixDouble > fieldGradPtr
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode runProgram()
boost::shared_ptr< EdgeEle > boundaryResidualVectorPipeline
boost::shared_ptr< FaceEle > domainResidualVectorPipeline
MoFEMErrorCode setupProblem()
boost::shared_ptr< DataAtGaussPoints > previousUpdate
boost::shared_ptr< FaceEle > domainTangentMatrixPipeline