v0.13.2
Loading...
Searching...
No Matches
mixed_poisson.cpp

MixedPoisson intended to show how to solve mixed formulation of the Dirichlet problem for the Poisson equation using error indicators and p-adaptivity

/**
* \file mixed_poisson.cpp
* \example mixed_poisson.cpp
*
* MixedPoisson intended to show how to solve mixed formulation of the Dirichlet
* problem for the Poisson equation using error indicators and p-adaptivity
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
GAUSS>::OpMixDivTimesScalar<2>;
inline double sqr(double x) { return x * x; }
struct MixedPoisson {
MixedPoisson(MoFEM::Interface &m_field) : mField(m_field) {}
private:
int baseOrder;
//! [Analytical function]
static double analyticalFunction(const double x, const double y,
const double z) {
return exp(-100. * (sqr(x) + sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
}
//! [Analytical function]
//! [Analytical function gradient]
static VectorDouble analyticalFunctionGrad(const double x, const double y,
const double z) {
res.resize(2);
res[0] = -exp(-100. * (sqr(x) + sqr(y))) *
(200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
res[1] = -exp(-100. * (sqr(x) + sqr(y))) *
(200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
return res;
}
//! [Analytical function gradient]
//! [Source function]
static double sourceFunction(const double x, const double y, const double z) {
return -exp(-100. * (sqr(x) + sqr(y))) *
(400. * M_PI *
(x * cos(M_PI * y) * sin(M_PI * x) +
y * cos(M_PI * x) * sin(M_PI * y)) +
2. * (20000. * (sqr(x) + sqr(y)) - 200. - sqr(M_PI)) *
cos(M_PI * x) * cos(M_PI * y));
}
//! [Source function]
const char *name, DataType type,
Tag &tag_handle) {
int int_val = 0;
double double_val = 0;
switch (type) {
case MB_TYPE_INTEGER:
CHKERR m_field.get_moab().tag_get_handle(
name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
break;
case MB_TYPE_DOUBLE:
CHKERR m_field.get_moab().tag_get_handle(
name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
break;
default:
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong data type %d for tag", type);
}
}
MoFEMErrorCode checkError(int iter_num = 0);
MoFEMErrorCode outputResults(int iter_num = 0);
struct CommonData {
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxValsGrad;
boost::shared_ptr<MatrixDouble> approxFlux;
enum VecElements {
};
};
boost::shared_ptr<CommonData> commonDataPtr;
struct OpError : public DomainEleOp {
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr,
MoFEM::Interface &m_field)
: DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr),
mField(m_field) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
};
};
//! [Run programme]
}
//! [Run programme]
//! [Read mesh]
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
//! [Read mesh]
//! [Set up problem]
// Note that in 2D case HDIV and HCURL spaces are isomorphic, and therefore
// only base for HCURL has been implemented in 2D. Base vectors for HDIV space
// are be obtained after rotation of HCURL base vectors by a right angle
1);
// We use AINSWORTH_LEGENDRE_BASE since DEMKOWICZ_JACOBI_BASE for triangle
// is not yet implemented for L2 space. For quads DEMKOWICZ_JACOBI_BASE and
// AINSWORTH_LEGENDRE_BASE are construcreed in the same way
CHKERR simpleInterface->addDomainField("U", L2, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ref_iter_num", &refIterNum,
PETSC_NULL);
baseOrder = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-base_order", &baseOrder,
PETSC_NULL);
CHKERR simpleInterface->setFieldOrder("FLUX", baseOrder);
CHKERR simpleInterface->setFieldOrder("U", baseOrder - 1);
CHKERR mField.get_moab().get_entities_by_dimension(0, 2, domainEntities,
false);
Tag th_order;
CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
for (auto ent : domainEntities) {
CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &baseOrder);
}
}
//! [Set up problem]
//! [Set integration rule]
auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
}
//! [Set integration rule]
//! [Create common data]
commonDataPtr = boost::make_shared<CommonData>();
PetscInt ghosts[4] = {0, 1, 2, 3};
commonDataPtr->petscVec =
createSmartGhostVector(mField.get_comm(), 4, 4, 0, ghosts);
else
commonDataPtr->petscVec =
createSmartGhostVector(mField.get_comm(), 0, 4, 4, ghosts);
commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
}
//! [Create common data]
//! [Assemble system]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainLhsPipeline().clear();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacHcurlFace(inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
auto beta = [](const double, const double, const double) { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHdivHdiv("FLUX", "FLUX", beta));
auto unity = []() { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHdivU("FLUX", "U", unity, true));
auto source = [&](const double x, const double y, const double z) {
return -sourceFunction(x, y, z);
};
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainSource("U", source));
}
//! [Assemble system]
//! [Solve]
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simpleInterface->getDM();
auto D = smartCreateDMVector(dm);
CHKERR VecZeroEntries(D);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
//! [Refine]
Tag th_error_ind, th_order;
CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE, th_error_ind);
CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
std::vector<Range> refinement_levels;
refinement_levels.resize(refIterNum + 1);
for (auto ent : domainEntities) {
double err_indic = 0;
CHKERR mField.get_moab().tag_get_data(th_error_ind, &ent, 1, &err_indic);
int order, new_order;
CHKERR mField.get_moab().tag_get_data(th_order, &ent, 1, &order);
new_order = order + 1;
Range refined_ents;
refined_ents.insert(ent);
Range adj;
CHKERR mField.get_moab().get_adjacencies(&ent, 1, 1, false, adj,
moab::Interface::UNION);
refined_ents.merge(adj);
refinement_levels[new_order - baseOrder].merge(refined_ents);
CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &new_order);
}
}
for (int ll = 1; ll < refinement_levels.size(); ll++) {
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
refinement_levels[ll]);
CHKERR mField.set_field_order(refinement_levels[ll], "FLUX",
baseOrder + ll);
CHKERR mField.set_field_order(refinement_levels[ll], "U",
baseOrder + ll - 1);
}
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("FLUX");
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
}
//! [Refine]
//! [Solve and refine loop]
for (int ii = 0; ii < refIterNum; ii++) {
CHKERR checkError(ii + 1);
}
}
//! [Solve and refine loop]
//! [Check error]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacHcurlFace(inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacL2ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
commonDataPtr->approxValsGrad));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHVecVectorField<3>("FLUX", commonDataPtr->approxFlux));
pipeline_mng->getOpDomainRhsPipeline().push_back(
CHKERR VecZeroEntries(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(commonDataPtr->petscVec, INSERT_VALUES,
SCATTER_FORWARD);
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Global error L2 norm: " << std::sqrt(array[0]);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Global error H1 seminorm: " << std::sqrt(array[1]);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Global error indicator: " << std::sqrt(array[2]);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Total number of elements: " << (int)array[3];
totalElementNumber = (int)array[3];
CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
std::vector<Tag> tag_handles;
tag_handles.resize(4);
CHKERR getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE, tag_handles[0]);
CHKERR getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
tag_handles[1]);
CHKERR getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
tag_handles[2]);
CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, tag_handles[3]);
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
if (pcomm == NULL)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Communicator not set");
tag_handles.push_back(pcomm->part_tag());
std::ostringstream strm;
strm << "error_" << iter_num << ".h5m";
CHKERR mField.get_moab().write_file(strm.str().c_str(), "MOAB",
"PARALLEL=WRITE_PART", 0, 0,
tag_handles.data(), tag_handles.size());
}
//! [Check error]
//! [Output results]
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe =
boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOJacForFace(jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacHcurlFace(inv_jac_ptr));
auto u_ptr = boost::make_shared<VectorDouble>();
auto flux_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", u_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHVecVectorField<3>("FLUX", flux_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"U", u_ptr}},
OpPPMap::DataMapMat{{"FLUX", flux_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
std::ostringstream strm;
strm << "out_" << iter_num << ".h5m";
CHKERR post_proc_fe->writeFile(strm.str().c_str());
}
//! [Output results]
//! [OpError]
EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
auto t_val_grad = getFTensor1FromMat<2>(*(commonDataPtr->approxValsGrad));
auto t_flux = getFTensor1FromMat<3>(*(commonDataPtr->approxFlux));
auto t_coords = getFTensor1CoordsAtGaussPts();
double error_l2 = 0;
double error_h1 = 0;
double error_ind = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * area;
double diff = t_val - MixedPoisson::analyticalFunction(
t_coords(0), t_coords(1), t_coords(2));
error_l2 += alpha * sqr(diff);
t_coords(0), t_coords(1), t_coords(2));
auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
t_diff(i) = t_val_grad(i) - t_fun_grad(i);
error_h1 += alpha * t_diff(i) * t_diff(i);
t_diff(i) = t_val_grad(i) - t_flux(i);
error_ind += alpha * t_diff(i) * t_diff(i);
++t_w;
++t_val;
++t_val_grad;
++t_flux;
++t_coords;
}
Tag th_error_l2, th_error_h1, th_error_ind;
CHKERR MixedPoisson::getTagHandle(mField, "ERROR_L2_NORM", MB_TYPE_DOUBLE,
th_error_l2);
CHKERR MixedPoisson::getTagHandle(mField, "ERROR_H1_SEMINORM", MB_TYPE_DOUBLE,
th_error_h1);
CHKERR MixedPoisson::getTagHandle(mField, "ERROR_INDICATOR", MB_TYPE_DOUBLE,
th_error_ind);
CHKERR mField.get_moab().tag_set_data(th_error_l2, &ent, 1, &error_l2);
CHKERR mField.get_moab().tag_set_data(th_error_h1, &ent, 1, &error_h1);
CHKERR mField.get_moab().tag_set_data(th_error_ind, &ent, 1, &error_ind);
constexpr std::array<int, 4> indices = {
std::array<double, 4> values;
values[0] = error_l2;
values[1] = error_h1;
values[2] = error_ind;
values[3] = 1.;
CHKERR VecSetValues(commonDataPtr->petscVec, 4, indices.data(), values.data(),
ADD_VALUES);
}
//! [OpError]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example problem
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "MixedPoisson");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [MixedPoisson]
MixedPoisson ex(m_field);
CHKERR ex.runProblem();
//! [MixedPoisson]
}
}
static Index< 'p', 3 > p
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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:511
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMoFEM.cpp:1269
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
double D
double sqr(double x)
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
boost::shared_ptr< MatrixDouble > approxFlux
SmartPetscObj< Vec > petscVec
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > approxValsGrad
MoFEM::Interface & mField
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Output results]
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode solveRefineLoop()
[Refine]
MoFEMErrorCode assembleSystem()
[Create common data]
Simple * simpleInterface
Range domainEntities
MoFEMErrorCode outputResults(int iter_num=0)
[Check error]
MoFEMErrorCode refineOrder()
[Solve]
MoFEMErrorCode createCommonData()
[Set integration rule]
MoFEMErrorCode checkError(int iter_num=0)
[Solve and refine loop]
static double sourceFunction(const double x, const double y, const double z)
[Analytical function gradient]
MoFEMErrorCode solveSystem()
[Assemble system]
static VectorDouble analyticalFunctionGrad(const double x, const double y, const double z)
[Analytical function]
MoFEM::Interface & mField
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode setIntegrationRules()
[Set up problem]
static double analyticalFunction(const double x, const double y, const double z)
[Analytical function]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setupProblem()
[Read mesh]
double errorIndicatorIntegral
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
MoFEMErrorCode readMesh()
[Run programme]
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() 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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
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
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get vector field for H-div approximation.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]