14static char help[] =
"...\n\n";
25 GAUSS>::OpMixDivTimesScalar<2>;
29inline double sqr(
double x) {
return x * x; }
51 return exp(-100. * (
sqr(x) +
sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
60 res[0] = -exp(-100. * (
sqr(x) +
sqr(y))) *
61 (200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
62 res[1] = -exp(-100. * (
sqr(x) +
sqr(y))) *
63 (200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
69 static double sourceFunction(
const double x,
const double y,
const double z) {
70 return -exp(-100. * (
sqr(x) +
sqr(y))) *
72 (x * cos(M_PI * y) * sin(M_PI * x) +
73 y * cos(M_PI * x) * sin(M_PI * y)) +
74 2. * (20000. * (
sqr(x) +
sqr(y)) - 200. -
sqr(M_PI)) *
75 cos(M_PI * x) * cos(M_PI * y));
80 const char *name, DataType type,
84 double double_val = 0;
88 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
92 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
96 "Wrong data type %d for tag", type);
132 OpError(boost::shared_ptr<CommonData> &common_data_ptr,
210 auto rule = [](int, int,
int p) ->
int {
return 2 *
p; };
224 PetscInt ghosts[3] = {0, 1, 2};
231 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
232 commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
233 commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
252 auto unity = []() {
return 1; };
254 new OpHdivU(
"Q",
"U", unity,
true));
255 auto source = [&](
const double x,
const double y,
const double z) {
269 CHKERR KSPSetFromOptions(solver);
278 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
279 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
288 Tag th_error_ind, th_order;
292 std::vector<Range> refinement_levels;
293 refinement_levels.resize(iter_num + 1);
295 double err_indic = 0;
298 int order, new_order;
300 new_order =
order + 1;
303 refined_ents.insert(ent);
306 moab::Interface::UNION);
307 refined_ents.merge(adj);
308 refinement_levels[new_order -
initOrder].merge(refined_ents);
313 for (
int ll = 1; ll < refinement_levels.size(); ll++) {
315 refinement_levels[ll]);
319 <<
"setting approximation order higher than 8 is not currently "
352 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Refinement iteration " << iter_num;
364 "Too many refinement iterations");
407 <<
"Global error indicator (max): " <<
commonDataPtr->maxErrorIndicator;
409 <<
"Global error indicator (total): "
412 <<
"Global error L2 norm: "
415 <<
"Global error H1 seminorm: "
421 std::vector<Tag> tag_handles;
422 tag_handles.resize(4);
430 ParallelComm *pcomm =
435 tag_handles.push_back(pcomm->part_tag());
436 std::ostringstream strm;
437 strm <<
"error_" << iter_num <<
".h5m";
439 "PARALLEL=WRITE_PART", 0, 0,
440 tag_handles.data(), tag_handles.size());
453 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
457 auto u_ptr = boost::make_shared<VectorDouble>();
458 auto flux_ptr = boost::make_shared<MatrixDouble>();
459 post_proc_fe->getOpPtrVector().push_back(
461 post_proc_fe->getOpPtrVector().push_back(
466 post_proc_fe->getOpPtrVector().push_back(
468 new OpPPMap(post_proc_fe->getPostProcMesh(),
469 post_proc_fe->getMapGaussPts(),
483 pipeline_mng->getDomainRhsFE() = post_proc_fe;
484 CHKERR pipeline_mng->loopFiniteElements();
486 std::ostringstream strm;
487 strm <<
"out_" << iter_num <<
".h5m";
488 CHKERR post_proc_fe->writeFile(strm.str().c_str());
497 const int nb_integration_pts =
getGaussPts().size2();
501 auto t_val_grad = getFTensor1FromMat<2>(*(
commonDataPtr->approxValsGrad));
502 auto t_flux = getFTensor1FromMat<3>(*(
commonDataPtr->approxFlux));
509 double error_ind = 0;
510 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
511 const double alpha = t_w * area;
514 t_coords(0), t_coords(1), t_coords(2));
515 error_l2 += alpha *
sqr(diff);
518 t_coords(0), t_coords(1), t_coords(2));
519 auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
520 t_diff(
i) = t_val_grad(
i) - t_fun_grad(
i);
521 error_h1 += alpha * t_diff(
i) * t_diff(
i);
523 t_diff(
i) = t_val_grad(
i) - t_flux(
i);
524 error_ind += alpha * t_diff(
i) * t_diff(
i);
534 Tag th_error_l2, th_error_h1, th_error_ind;
542 double error_l2_norm = sqrt(error_l2);
543 double error_h1_seminorm = sqrt(error_h1);
544 double error_ind_local = sqrt(error_ind);
555 constexpr std::array<int, CommonData::LAST_ELEMENT> indices = {
558 std::array<double, CommonData::LAST_ELEMENT> values;
559 values[0] = error_l2;
560 values[1] = error_h1;
561 values[2] = error_ind;
563 indices.data(), values.data(), ADD_VALUES);
568int main(
int argc,
char *argv[]) {
574 auto core_log = logging::core::get();
582 DMType dm_name =
"DMMOFEM";
587 moab::Core mb_instance;
588 moab::Interface &moab = mb_instance;
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMSetUp_MoFEM(DM dm)
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.
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 2 > OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
auto createGhostVector(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]
OpError(boost::shared_ptr< CommonData > &common_data_ptr, MoFEM::Interface &m_field)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode solveRefineLoop()
[Refine]
MoFEMErrorCode assembleSystem()
[Create common data]
MoFEMErrorCode outputResults(int iter_num=0)
[Check error]
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
MixedPoisson(MoFEM::Interface &m_field)
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 refineOrder(int iter_num=0)
[Solve]
MoFEMErrorCode readMesh()
[Run programme]
Add operators pushing bases from local to physical configuration.
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
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.
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
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.
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.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
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]