v0.14.0
SCL-3: Poisson's equation (Lagrange multiplier)
Note
Prerequisite of this tutorial is SCL-2: Poisson's equation (non-homogeneous BC)


Note
Intended learning outcome:
  • an alternative way to handle non-homogeneous boundary condition using an additional field of Lagrange multipliers
  • idea of PETSc PCFIELDSPLIT block solver and how to implement it
  • adding the additional field to the output file (postprocessing)

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 __POISSON2DLAGRANGEMULTIPLIER_HPP__
#define __POISSON2DLAGRANGEMULTIPLIER_HPP__
#include <stdlib.h>
// const double body_source = 10;
typedef boost::function<double(const double, const double, const double)>
struct OpDomainLhsK : public OpFaceEle {
public:
OpDomainLhsK(std::string row_field_name, std::string col_field_name,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
: OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
boundaryMarker(boundary_marker) {
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();
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 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 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) += t_row_diff_base(i) * t_col_diff_base(i) * 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;
}
// FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
ADD_VALUES);
// Fill values of symmetric stiffness matrix in global system of equations
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocLhs(0, 0),
ADD_VALUES);
}
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
struct OpDomainRhsF : public OpFaceEle {
public:
OpDomainRhsF(std::string field_name, ScalarFunc source_term_function,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
sourceTermFunc(source_term_function), 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 base functions
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 * area;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locRhs[rr] += t_base * body_source * 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;
}
// FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
CHKERR VecSetValues(getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
struct OpBoundaryLhsC : public OpEdgeEle {
public:
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
: OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
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) {
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 functions 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 functions 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 functions on column
++t_col_base;
}
// move to the next base functions 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(getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0),
ADD_VALUES);
}
}
private:
};
struct OpBoundaryRhsG : public OpEdgeEle {
public:
OpBoundaryRhsG(std::string field_name, ScalarFunc boundary_function)
boundaryFunc(boundary_function) {}
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 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));
for (int rr = 0; rr != nb_dofs; rr++) {
locBoundaryRhs[rr] += t_base * boundary_term * 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;
}
// FILL VALUES OF LOCAL VECTOR ENTRIES TO THE GLOBAL VECTOR
ADD_VALUES);
}
}
private:
};
}; // namespace Poisson2DLagrangeMultiplierOperators
#endif //__POISSON2DLAGRANGEMULTIPLIER_HPP__

Implementation of the main program (*.cpp)

#include <stdlib.h>
using namespace MoFEM;
static char help[] = "...\n\n";
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 boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
// 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
Simple *simpleInterface;
// Field name and approximation order
std::string domainField; // displacement field
std::string boundaryField; // Lagrange multiplier field
int oRder;
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
// Boundary entities marked for fieldsplit (block) solver - optional
Range boundaryEntitiesForFieldsplit;
};
MoFEM::Interface &m_field)
: domainField("U"), boundaryField("L"), mField(m_field) {}
}
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
// Get boundary edges marked in block named "BOUNDARY_CONDITION"
auto get_ents_on_mesh = [&]() {
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());
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
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 LHS of
// domain elements
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJac<2>(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// domain elements
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating LHS of
// boundary elements (Lagrange multiplier)
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating RHS of
// boundary elements (Lagrange multiplier)
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
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); };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
// Create RHS and solution vectors
auto dm = simpleInterface->getDM();
auto F = createDMVector(dm);
auto D = vectorDuplicate(F);
// Setup fieldsplit (block) solver - optional: yes/no
if (1) {
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
// Set up FIELDSPLIT, only when user set -pc_type fieldsplit
// Identify the index for boundary entities, remaining will be for domain
// Then split the fields for boundary and domain for solving
if (is_pcfs == PETSC_TRUE) {
IS is_domain, is_boundary;
cerr << "Running FIELDSPLIT..." << endl;
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 KSPSetUp(ksp_solver);
// Solve the system
CHKERR KSPSolve(ksp_solver, F, D);
// Scatter result data on the mesh
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto d_ptr = boost::make_shared<VectorDouble>();
auto l_ptr = boost::make_shared<VectorDouble>();
auto post_proc_domain_fe = boost::make_shared<PostProcFaceEle>(mField);
post_proc_domain_fe->getOpPtrVector().push_back(
post_proc_domain_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_domain_fe->getPostProcMesh(),
post_proc_domain_fe->getMapGaussPts(), {{domainField, d_ptr}},
{}, {}, {}));
pipeline_mng->getDomainRhsFE() = post_proc_domain_fe;
auto post_proc_boundary_fe = boost::make_shared<PostProcEdgeEle>(mField);
post_proc_boundary_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(boundaryField, l_ptr));
post_proc_boundary_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_boundary_fe->getPostProcMesh(),
post_proc_boundary_fe->getMapGaussPts(),
{{boundaryField, l_ptr}}, {}, {}, {}));
pipeline_mng->getBoundaryRhsFE() = post_proc_boundary_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_domain_fe->writeFile("out_result_domain.h5m");
CHKERR post_proc_boundary_fe->writeFile("out_result_boundary.h5m");
}
}
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);
// 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
Poisson2DLagrangeMultiplier poisson_problem(m_field);
CHKERR poisson_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPf
Vec getKSPf() const
Definition: ForcesAndSourcesCore.hpp:1087
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
Poisson2DLagrangeMultiplier::sourceTermFunction
static double sourceTermFunction(const double x, const double y, const double z)
Definition: poisson_2d_lagrange_multiplier.cpp:33
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF::OpDomainRhsF
OpDomainRhsF(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
Definition: poisson_2d_lagrange_multiplier.hpp:106
MoFEM::PostProcBrokenMeshInMoab< EdgeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:684
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
Poisson2DLagrangeMultiplier::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: poisson_2d_lagrange_multiplier.cpp:134
Poisson2DLagrangeMultiplier::boundaryField
std::string boundaryField
Definition: poisson_2d_lagrange_multiplier.cpp:51
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC::doWork
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.
Definition: poisson_2d_lagrange_multiplier.hpp:176
Poisson2DLagrangeMultiplier::solveSystem
MoFEMErrorCode solveSystem()
Definition: poisson_2d_lagrange_multiplier.cpp:197
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC::transLocBoundaryLhs
MatrixDouble transLocBoundaryLhs
Definition: poisson_2d_lagrange_multiplier.hpp:233
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
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::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1239
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1274
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1235
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK::transLocLhs
MatrixDouble transLocLhs
Definition: poisson_2d_lagrange_multiplier.hpp:101
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG::OpBoundaryRhsG
OpBoundaryRhsG(std::string field_name, ScalarFunc boundary_function)
Definition: poisson_2d_lagrange_multiplier.hpp:238
Poisson2DLagrangeMultiplier::Poisson2DLagrangeMultiplier
Poisson2DLagrangeMultiplier(MoFEM::Interface &m_field)
Definition: poisson_2d_lagrange_multiplier.cpp:61
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
Poisson2DLagrangeMultiplierOperators::i
FTensor::Index< 'i', 2 > i
Definition: poisson_2d_lagrange_multiplier.hpp:17
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
Poisson2DLagrangeMultiplier::readMesh
MoFEMErrorCode readMesh()
Definition: poisson_2d_lagrange_multiplier.cpp:65
Poisson2DLagrangeMultiplier::outputResults
MoFEMErrorCode outputResults()
Definition: poisson_2d_lagrange_multiplier.cpp:249
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK::OpDomainLhsK
OpDomainLhsK(std::string row_field_name, std::string col_field_name, boost::shared_ptr< std::vector< unsigned char >> boundary_marker=nullptr)
Definition: poisson_2d_lagrange_multiplier.hpp:25
Poisson2DLagrangeMultiplier::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: poisson_2d_lagrange_multiplier.cpp:179
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK::locLhs
MatrixDouble locLhs
Definition: poisson_2d_lagrange_multiplier.hpp:101
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_lagrange_multiplier.hpp:165
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
double
convert.type
type
Definition: convert.py:64
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC
Definition: poisson_2d_lagrange_multiplier.hpp:169
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Poisson2DLagrangeMultiplier::boundaryEntitiesForFieldsplit
Range boundaryEntitiesForFieldsplit
Definition: poisson_2d_lagrange_multiplier.cpp:58
Poisson2DHomogeneousOperators::body_source
const double body_source
Definition: poisson_2d_homogeneous.hpp:36
Poisson2DLagrangeMultiplier::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: poisson_2d_lagrange_multiplier.cpp:39
Poisson2DLagrangeMultiplier::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: poisson_2d_lagrange_multiplier.cpp:93
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:1201
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
Poisson2DLagrangeMultiplier::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_lagrange_multiplier.cpp:47
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: poisson_2d_lagrange_multiplier.hpp:242
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
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:3200
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG::locBoundaryRhs
VectorDouble locBoundaryRhs
Definition: poisson_2d_lagrange_multiplier.hpp:296
Poisson2DLagrangeMultiplier::setupProblem
MoFEMErrorCode setupProblem()
Definition: poisson_2d_lagrange_multiplier.cpp:75
Poisson2DLagrangeMultiplier::domainField
std::string domainField
Definition: poisson_2d_lagrange_multiplier.cpp:50
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_lagrange_multiplier.hpp:100
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF
Definition: poisson_2d_lagrange_multiplier.hpp:104
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 >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC::locBoundaryLhs
MatrixDouble locBoundaryLhs
Definition: poisson_2d_lagrange_multiplier.hpp:233
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF::sourceTermFunc
ScalarFunc sourceTermFunc
Definition: poisson_2d_lagrange_multiplier.hpp:164
Poisson2DLagrangeMultiplierOperators
Definition: poisson_2d_lagrange_multiplier.hpp:15
poisson_2d_lagrange_multiplier.hpp
Range
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::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::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:372
Poisson2DLagrangeMultiplier
Definition: poisson_2d_lagrange_multiplier.cpp:15
_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:1102
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG::boundaryFunc
ScalarFunc boundaryFunc
Definition: poisson_2d_lagrange_multiplier.hpp:295
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: poisson_2d_lagrange_multiplier.hpp:111
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
Poisson2DLagrangeMultiplier::runProgram
MoFEMErrorCode runProgram()
Definition: poisson_2d_lagrange_multiplier.cpp:286
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
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:282
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Poisson2DLagrangeMultiplierOperators::ScalarFunc
boost::function< double(const double, const double, const double)> ScalarFunc
Definition: poisson_2d_lagrange_multiplier.hpp:21
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
Poisson2DLagrangeMultiplier::oRder
int oRder
Definition: poisson_2d_lagrange_multiplier.cpp:52
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK::doWork
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.
Definition: poisson_2d_lagrange_multiplier.hpp:32
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
Poisson2DLagrangeMultiplier::mField
MoFEM::Interface & mField
Definition: poisson_2d_lagrange_multiplier.cpp:46
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF::locRhs
VectorDouble locRhs
Definition: poisson_2d_lagrange_multiplier.hpp:166
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
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:416
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC::OpBoundaryLhsC
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
Definition: poisson_2d_lagrange_multiplier.hpp:171
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1103
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG
Definition: poisson_2d_lagrange_multiplier.hpp:236
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK
Definition: poisson_2d_lagrange_multiplier.hpp:23
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
Definition: ForcesAndSourcesCore.hpp:1268
Poisson2DLagrangeMultiplier::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_lagrange_multiplier.cpp:55
F
@ F
Definition: free_surface.cpp:394
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698