v0.14.0
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>
#include <MoFEM.hpp>
using namespace MoFEM;
PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
typedef boost::function<double(const double, const double, const double)>
struct OpDomainLhs : public AssemblyDomainEleOp {
public:
OpDomainLhs(
std::string row_field_name, std::string col_field_name,
boost::shared_ptr<VectorDouble> field_vec,
boost::shared_ptr<MatrixDouble> field_grad_mat)
: AssemblyDomainEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
fieldVec(field_vec), fieldGradMat(field_grad_mat) {}
MoFEMErrorCode iNtegrate(EntData &row_data,
EntData &col_data) {
const int nb_row_dofs = row_data.getIndices().size();
const int nb_col_dofs = col_data.getIndices().size();
// 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 solution (field value) at integration points
auto t_field = getFTensor0FromVec(*fieldVec);
// get gradient of field at integration points
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
// 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;
}
}
private:
boost::shared_ptr<VectorDouble> fieldVec;
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
struct OpDomainRhs : public AssemblyDomainEleOp {
public:
OpDomainRhs(
std::string field_name, ScalarFunc source_term_function,
boost::shared_ptr<VectorDouble> field_vec,
boost::shared_ptr<MatrixDouble> field_grad_mat)
sourceTermFunc(source_term_function), fieldVec(field_vec),
fieldGradMat(field_grad_mat) {}
MoFEMErrorCode iNtegrate(EntData &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 coordinates of the integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// get solution (field value) at integration point
auto t_field = getFTensor0FromVec(*fieldVec);
// get gradient of field value of integration point
auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
// 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 != AssemblyDomainEleOp::nbRows; rr++) {
nf[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;
}
}
private:
ScalarFunc sourceTermFunc;
boost::shared_ptr<VectorDouble> fieldVec;
boost::shared_ptr<MatrixDouble> fieldGradMat;
};
}; // namespace NonlinearPoissonOps
#endif //__NONLINEARPOISSON2D_HPP__

Implementation of the main program (*.cpp)

#include <stdlib.h>
#include <MoFEM.hpp>
using namespace MoFEM;
using namespace NonlinearPoissonOps;
constexpr int SPACE_DIM = 2;
static char help[] = "...\n\n";
inline double sqr(const double x) { return x * x; }
inline double cube(const double x) { return x * x * x; }
public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
// Function to calculate the Source term
static double sourceTermFunction(const double x, const double y,
const double z) {
return 2 * M_PI * M_PI *
(cos(M_PI * x) * cos(M_PI * y) +
cube(cos(M_PI * x)) * cube(cos(M_PI * y)) -
cos(M_PI * x) * cos(M_PI * y) *
(sqr(sin(M_PI * x)) * sqr(cos(M_PI * y)) +
sqr(sin(M_PI * y)) * sqr(cos(M_PI * x))));
}
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y,
const double z) {
return -cos(M_PI * x) *
cos(M_PI * y); // here should put the negative of the proper formula
}
// Main interfaces
Simple *simpleInterface;
// Field name and approximation order
std::string domainField = "POTENTIAL";
int order;
// PETSc vector for storing norms
SmartPetscObj<Vec> errorVec, exactVec;
int atomTest = 0;
enum NORMS { NORM = 0, LAST_NORM };
enum EX_NORMS { EX_NORM = 0, LAST_EX_NORM };
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
// Object needed for postprocessing
boost::shared_ptr<PostProcEle> postProc;
// Boundary entities marked for fieldsplit (block) solver - optional
Range boundaryEntitiesForFieldsplit;
};
: mField(m_field) {}
}
}
Range domain_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
true);
auto get_ents_by_dim = [&](const auto dim) {
if (dim == SPACE_DIM) {
return domain_ents;
} else {
Range ents;
if (dim == 0)
CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
else
CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
return ents;
}
};
// Select base
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(SPACE_DIM);
if (domain_ents.empty())
const auto type = type_from_handle(domain_ents[0]);
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
}
return NOBASE;
};
auto base = get_base();
order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atomTest,
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; };
auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
auto pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_rhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(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>>();
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 pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR AddHOOps<2, 2, 2>::add(pipeline_mng->getOpDomainLhsPipeline(), {H1});
CHKERR AddHOOps<2, 2, 2>::add(pipeline_mng->getOpDomainRhsPipeline(), {H1});
auto add_domain_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc(domainField, true, boundaryMarker));
auto data_u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldValues(domainField, data_u_at_gauss_pts));
pipeline.push_back(new OpCalculateScalarFieldGradient<2>(
domainField, grad_u_at_gauss_pts));
pipeline.push_back(new OpDomainLhs(
domainField, domainField, data_u_at_gauss_pts, grad_u_at_gauss_pts));
pipeline.push_back(new OpUnSetBc(domainField));
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc(domainField, true, boundaryMarker));
auto data_u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateScalarFieldValues(domainField, data_u_at_gauss_pts));
pipeline.push_back(new OpCalculateScalarFieldGradient<2>(
domainField, grad_u_at_gauss_pts));
pipeline.push_back(new OpDomainRhs(domainField, sourceTermFunction,
data_u_at_gauss_pts,
grad_u_at_gauss_pts));
pipeline.push_back(new OpUnSetBc(domainField));
};
auto add_boundary_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc(domainField, false, boundaryMarker));
pipeline.push_back(new OpBoundaryLhs(
[](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc(domainField));
};
auto add_boundary_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc(domainField, false, boundaryMarker));
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(
pipeline.push_back(new OpBoundaryRhs(
domainField, u_at_gauss_pts,
[](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc(domainField));
};
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_rhs_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);
// Set up FIELDSPLIT, only when option used -pc_type fieldsplit
if (is_pcfs == PETSC_TRUE) {
auto name_prb = simple->getProblemName();
SmartPetscObj<IS> is_all_bc;
CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
name_prb, ROW, domainField, 0, 1, is_all_bc,
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 = vectorDuplicate(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(
post_proc->getOpPtrVector().push_back(
new OpPPMap(post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
{{domainField, 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");
}
//! [Check]
auto check_result_fe_ptr = boost::make_shared<DomainEle>(mField);
auto errorVec =
check_result_fe_ptr->getOpPtrVector(), {H1})),
"Apply transform");
check_result_fe_ptr->getRuleHook = [](int, int, int p) { return p; };
auto analyticalFunction = [&](double x, double y, double z) {
return cos(M_PI * x) * cos(M_PI * y);
};
auto u_ptr = boost::make_shared<VectorDouble>();
check_result_fe_ptr->getOpPtrVector().push_back(
auto mValFuncPtr = boost::make_shared<VectorDouble>();
check_result_fe_ptr->getOpPtrVector().push_back(
new OpGetTensor0fromFunc(mValFuncPtr, analyticalFunction));
check_result_fe_ptr->getOpPtrVector().push_back(
new OpCalcNormL2Tensor0(u_ptr, errorVec, NORM, mValFuncPtr));
check_result_fe_ptr->getOpPtrVector().push_back(
CHKERR VecZeroEntries(errorVec);
CHKERR VecZeroEntries(exactVec);
check_result_fe_ptr);
CHKERR VecAssemblyBegin(errorVec);
CHKERR VecAssemblyEnd(errorVec);
CHKERR VecAssemblyBegin(exactVec);
CHKERR VecAssemblyEnd(exactVec);
MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
// print norm in general log
if (mField.get_comm_rank() == 0) {
const double *L2_norms;
const double *Ex_norms;
CHKERR VecGetArrayRead(errorVec, &L2_norms);
CHKERR VecGetArrayRead(exactVec, &Ex_norms);
MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
<< "L2_NORM: " << std::sqrt(L2_norms[NORM]);
MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
<< "NORMALISED_ERROR: "
<< (std::sqrt(L2_norms[NORM]) / std::sqrt(Ex_norms[EX_NORM]));
CHKERR VecRestoreArrayRead(errorVec, &L2_norms);
CHKERR VecRestoreArrayRead(exactVec, &Ex_norms);
}
// compare norm for ctest
const double *L2_norms;
const double *Ex_norms;
CHKERR VecGetArrayRead(errorVec, &L2_norms);
CHKERR VecGetArrayRead(exactVec, &Ex_norms);
double ref_percentage = 4.4e-04;
double normalised_error;
switch (atomTest) {
case 1: // 2D
normalised_error = std::sqrt(L2_norms[0]) / std::sqrt(Ex_norms[0]);
break;
default:
SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"atom test %d does not exist", atomTest);
}
if (normalised_error > ref_percentage) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test %d failed! Calculated Norm %3.16e is greater than "
"reference Norm %3.16e",
atomTest, normalised_error, ref_percentage);
}
CHKERR VecRestoreArrayRead(errorVec, &L2_norms);
CHKERR VecRestoreArrayRead(exactVec, &Ex_norms);
}
}
//! [Check]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
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
NonlinearPoisson poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
NonlinearPoisson::outputResults
MoFEMErrorCode outputResults()
Definition: nonlinear_poisson_2d.cpp:336
H1
@ H1
continuous field
Definition: definitions.h:85
NonlinearPoisson::EX_NORM
@ EX_NORM
Definition: nonlinear_poisson_2d.cpp:64
NonlinearPoisson::NonlinearPoisson
NonlinearPoisson(MoFEM::Interface &m_field)
Definition: nonlinear_poisson_2d.cpp:78
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NonlinearPoisson::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: nonlinear_poisson_2d.cpp:46
OpBoundaryRhs
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryRhs
Definition: nonlinear_poisson_2d.hpp:18
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
NOBASE
@ NOBASE
Definition: definitions.h:59
Poisson2DHomogeneousOperators::ScalarFunc
boost::function< double(const double, const double, const double)> ScalarFunc
Definition: poisson_2d_homogeneous.hpp:38
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
NonlinearPoisson::LAST_NORM
@ LAST_NORM
Definition: nonlinear_poisson_2d.cpp:63
NonlinearPoisson::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: nonlinear_poisson_2d.cpp:178
NonlinearPoisson::sourceTermFunction
static double sourceTermFunction(const double x, const double y, const double z)
Definition: nonlinear_poisson_2d.cpp:34
NonlinearPoisson::order
int order
Definition: nonlinear_poisson_2d.cpp:58
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
domainField
constexpr auto domainField
Definition: electrostatics.hpp:7
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
NonlinearPoissonOps
Definition: nonlinear_poisson_2d.hpp:31
NonlinearPoisson
Definition: nonlinear_poisson_2d.cpp:15
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
nonlinear_poisson_2d.hpp
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
NonlinearPoisson::NORM
@ NORM
Definition: nonlinear_poisson_2d.cpp:63
NonlinearPoisson::domainField
std::string domainField
Definition: nonlinear_poisson_2d.cpp:57
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::PipelineManager::EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: PipelineManager.hpp:36
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:251
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
NonlinearPoisson::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: nonlinear_poisson_2d.cpp:160
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
NonlinearPoisson::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: nonlinear_poisson_2d.cpp:67
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
NonlinearPoisson::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: nonlinear_poisson_2d.cpp:219
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
NonlinearPoisson::simpleInterface
Simple * simpleInterface
Definition: nonlinear_poisson_2d.cpp:54
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
BoundaryEleOp
OpBoundaryRhsSource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundaryRhsSource
Definition: nonlinear_poisson_2d.hpp:20
NonlinearPoisson::solveSystem
MoFEMErrorCode solveSystem()
Definition: nonlinear_poisson_2d.cpp:282
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
double
NonlinearPoisson::atomTest
int atomTest
Definition: nonlinear_poisson_2d.cpp:62
convert.type
type
Definition: convert.py:64
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
NonlinearPoisson::checkResults
MoFEMErrorCode checkResults()
[Check]
Definition: nonlinear_poisson_2d.cpp:369
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
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
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
NonlinearPoisson::LAST_EX_NORM
@ LAST_EX_NORM
Definition: nonlinear_poisson_2d.cpp:64
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
OpBoundaryLhs
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryLhs
Definition: nonlinear_poisson_2d.hpp:16
MoFEM::ProblemsManager::markDofs
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.
Definition: ProblemsManager.cpp:3548
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1898
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::OpGetTensor0fromFunc
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
Definition: NormsOperators.hpp:105
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
NonlinearPoisson::runProgram
MoFEMErrorCode runProgram()
Definition: nonlinear_poisson_2d.cpp:81
NonlinearPoisson::mField
MoFEM::Interface & mField
Definition: nonlinear_poisson_2d.cpp:53
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 2 >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
NonlinearPoisson::exactVec
SmartPetscObj< Vec > exactVec
Definition: nonlinear_poisson_2d.cpp:61
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
Range
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:236
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
NonlinearPoisson::errorVec
SmartPetscObj< Vec > errorVec
Definition: nonlinear_poisson_2d.cpp:61
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
MoFEM::OpCalcNormL2Tensor0
Get norm of input VectorDouble for Tensor0.
Definition: NormsOperators.hpp:15
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::Simple::addBoundaryField
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:354
sqr
double sqr(double x)
Definition: mixed_poisson.cpp:29
NonlinearPoisson::setupProblem
MoFEMErrorCode setupProblem()
Definition: nonlinear_poisson_2d.cpp:106
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::SmartPetscObj< Vec >
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:389
MoFEM::DMoFEMLoopFiniteElements
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:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
cube
double cube(const double x)
Definition: nonlinear_poisson_2d.cpp:13
NonlinearPoisson::readMesh
MoFEMErrorCode readMesh()
Definition: nonlinear_poisson_2d.cpp:96
NonlinearPoisson::boundaryEntitiesForFieldsplit
Range boundaryEntitiesForFieldsplit
Definition: nonlinear_poisson_2d.cpp:75
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698