static char help[] =
"...\n\n";
#include <BasicFiniteElements.hpp>
GAUSS>::OpMixDivTimesScalar<2>;
inline double sqr(
double x) {
return x * x; }
private:
Simple *simpleInterface;
Range domainEntities;
double errorIndicatorIntegral;
int totalElementNumber;
int baseOrder;
int refIterNum;
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);
}
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;
}
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));
}
const char *name, DataType
type,
Tag &tag_handle) {
int int_val = 0;
double double_val = 0;
case MB_TYPE_INTEGER:
name, 1,
type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
break;
case MB_TYPE_DOUBLE:
name, 1,
type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
break;
default:
"Wrong data type %d for tag",
type);
}
}
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxValsGrad;
boost::shared_ptr<MatrixDouble> approxFlux;
SmartPetscObj<Vec> petscVec;
enum VecElements {
ERROR_L2_NORM = 0,
ERROR_H1_SEMINORM,
ERROR_INDICATOR,
TOTAL_NUMBER,
LAST_ELEMENT
};
};
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr,
:
DomainEleOp(
"U", OPROW), commonDataPtr(common_data_ptr),
mField(m_field) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
};
};
}
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
1);
refIterNum = 0;
PETSC_NULL);
baseOrder = 2;
PETSC_NULL);
CHKERR simpleInterface->setFieldOrder(
"FLUX", baseOrder);
CHKERR simpleInterface->setFieldOrder(
"U", baseOrder - 1);
CHKERR simpleInterface->setUp();
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) {
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
PipelineManager *pipeline_mng = mField.
getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
}
commonDataPtr = boost::make_shared<CommonData>();
PetscInt ghosts[4] = {0, 1, 2, 3};
commonDataPtr->petscVec =
else
commonDataPtr->petscVec =
commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
}
PipelineManager *pipeline_mng = mField.
getInterface<PipelineManager>();
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(
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(
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
auto beta = [](const double, const double, const double) { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
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(
}
PipelineManager *pipeline_mng = mField.
getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
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);
Range refined_ents;
if (err_indic > errorIndicatorIntegral / totalElementNumber) {
refined_ents.insert(ent);
Range adj;
moab::Interface::UNION);
refined_ents.merge(adj);
refinement_levels[new_order - baseOrder].merge(refined_ents);
}
}
for (int ll = 1; ll < refinement_levels.size(); ll++) {
refinement_levels[ll]);
baseOrder + ll);
baseOrder + ll - 1);
}
mField.
getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
}
for (int ii = 0; ii < refIterNum; ii++) {
}
}
PipelineManager *pipeline_mng = mField.
getInterface<PipelineManager>();
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(
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(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacL2ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<2>("U",
commonDataPtr->approxValsGrad));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHVecVectorField<3>("FLUX", commonDataPtr->approxFlux));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpError(commonDataPtr, mField));
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);
<< "Global error L2 norm: " << std::sqrt(array[0]);
<< "Global error H1 seminorm: " << std::sqrt(array[1]);
<< "Global error indicator: " << std::sqrt(array[2]);
<<
"Total number of elements: " << (
int)array[3];
errorIndicatorIntegral = array[2];
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 =
if (pcomm == NULL)
tag_handles.push_back(pcomm->part_tag());
std::ostringstream strm;
strm << "error_" << iter_num << ".h5m";
"PARALLEL=WRITE_PART", 0, 0,
tag_handles.data(), tag_handles.size());
}
PipelineManager *pipeline_mng = mField.
getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe =
boost::make_shared<PostProcFaceOnRefinedMeshFor2D>(mField);
post_proc_fe->generateReferenceElementMesh();
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(
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(
post_proc_fe->addFieldValuesPostProc("FLUX");
post_proc_fe->addFieldValuesPostProc("U");
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());
}
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
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;
t_coords(0), t_coords(1), t_coords(2));
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;
th_error_l2);
th_error_h1);
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);
int index = CommonData::ERROR_L2_NORM;
constexpr std::array<int, 4> indices = {
CommonData::ERROR_L2_NORM, CommonData::ERROR_H1_SEMINORM,
CommonData::ERROR_INDICATOR, CommonData::TOTAL_NUMBER};
std::array<double, 4> values;
values[0] = error_l2;
values[1] = error_h1;
values[2] = error_ind;
values[3] = 1.;
ADD_VALUES);
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMSetUp_MoFEM(DM dm)
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.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
int main(int argc, char *argv[])
[OpError]
EntitiesFieldData::EntData EntData
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.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
auto createSmartGhostVector
Create smart ghost vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Output results]
MoFEMErrorCode solveRefineLoop()
[Refine]
MoFEMErrorCode assembleSystem()
[Create common data]
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]
MoFEMErrorCode solveSystem()
[Assemble system]
static VectorDouble analyticalFunctionGrad(const double x, const double y, const double z)
[Analytical function]
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]
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
MoFEMErrorCode readMesh()
[Run programme]
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)