#ifndef __POISSON2DLAGRANGEMULTIPLIER_HPP__
#define __POISSON2DLAGRANGEMULTIPLIER_HPP__
#include <stdlib.h>
using EntData = EntitiesFieldData::EntData;
typedef boost::function<
double(
const double,
const double,
const double)>
public:
OpDomainLhsK(std::string row_field_name, std::string col_field_name,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
}
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locLhs.resize(nb_row_dofs, nb_col_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += t_row_diff_base(
i) * t_col_diff_base(
i) *
a;
++t_col_diff_base;
}
++t_row_diff_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
ADD_VALUES);
}
}
}
private:
};
public:
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
if (nb_dofs) {
locRhs.resize(nb_dofs,
false);
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_dofs; rr++) {
++t_base;
}
++t_w;
++t_coords;
}
}
}
private:
};
public:
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
}
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
++t_col_base;
}
++t_row_base;
}
++t_w;
}
ADD_VALUES);
}
}
private:
};
public:
if (nb_dofs) {
const int nb_integration_points =
getGaussPts().size2();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
double boundary_term =
for (int rr = 0; rr != nb_dofs; rr++) {
++t_base;
}
++t_w;
++t_coords;
}
ADD_VALUES);
}
}
private:
};
};
#endif
#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.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
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.
boost::function< double(const double, const double, const double)> ScalarFunc
FTensor::Index< 'i', 2 > i
constexpr auto field_name
bool sYmm
If true assume that matrix is symmetric structure.
default operator for EDGE element
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 & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
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
MatrixDouble locBoundaryLhs
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.
MatrixDouble transLocBoundaryLhs
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
VectorDouble locBoundaryRhs
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpBoundaryRhsG(std::string field_name, ScalarFunc boundary_function)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpDomainRhsF(std::string field_name, ScalarFunc source_term_function, boost::shared_ptr< std::vector< unsigned char > > boundary_marker=nullptr)
ScalarFunc sourceTermFunc
#include <stdlib.h>
static char help[] =
"...\n\n";
public:
private:
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
}
const double z) {
return sin(x * 10.) * cos(y * 10.);
}
};
: domainField("U"), boundaryField("L"), mField(m_field) {}
}
}
auto get_ents_on_mesh = [&]() {
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
boundary_entities, true);
}
}
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](
Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
{
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
}
{
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
}
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);
}
CHKERR KSPSetFromOptions(ksp_solver);
if (1) {
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
IS is_domain, is_boundary;
cerr << "Running FIELDSPLIT..." << endl;
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
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(
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[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR poisson_problem.runProgram();
}
return 0;
}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#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
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
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Section manager is used to create indexes and sections.
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.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
std::string boundaryField
MoFEMErrorCode setupProblem()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode solveSystem()
static double sourceTermFunction(const double x, const double y, const double z)
Poisson2DLagrangeMultiplier(MoFEM::Interface &m_field)
MoFEMErrorCode boundaryCondition()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEM::Interface & mField
Range boundaryEntitiesForFieldsplit
MoFEMErrorCode outputResults()
static double boundaryFunction(const double x, const double y, const double z)
MoFEMErrorCode runProgram()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()