14static char help[] =
"...\n\n";
25 GAUSS>::OpMixDivTimesScalar<2>;
29inline double sqr(
double x) {
return x * x; }
49 return exp(-100. * (
sqr(x) +
sqr(y))) * cos(M_PI * x) * cos(M_PI * y);
58 res[0] = -exp(-100. * (
sqr(x) +
sqr(y))) *
59 (200. * x * cos(M_PI * x) + M_PI * sin(M_PI * x)) * cos(M_PI * y);
60 res[1] = -exp(-100. * (
sqr(x) +
sqr(y))) *
61 (200. * y * cos(M_PI * y) + M_PI * sin(M_PI * y)) * cos(M_PI * x);
67 static double sourceFunction(
const double x,
const double y,
const double z) {
68 return -exp(-100. * (
sqr(x) +
sqr(y))) *
70 (x * cos(M_PI * y) * sin(M_PI * x) +
71 y * cos(M_PI * x) * sin(M_PI * y)) +
72 2. * (20000. * (
sqr(x) +
sqr(y)) - 200. -
sqr(M_PI)) *
73 cos(M_PI * x) * cos(M_PI * y));
78 const char *name, DataType type,
82 double double_val = 0;
86 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
90 name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
94 "Wrong data type %d for tag", type);
129 OpError(boost::shared_ptr<CommonData> &common_data_ptr,
200 auto rule = [](int, int,
int p) ->
int {
return 2 *
p + 1; };
214 PetscInt ghosts[4] = {0, 1, 2, 3};
221 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
222 commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
223 commonDataPtr->approxFlux = boost::make_shared<MatrixDouble>();
237 auto det_ptr = boost::make_shared<VectorDouble>();
238 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
239 auto jac_ptr = boost::make_shared<MatrixDouble>();
254 auto unity = []() {
return 1; };
256 new OpHdivU(
"FLUX",
"U", unity,
true));
257 auto source = [&](
const double x,
const double y,
const double z) {
271 CHKERR KSPSetFromOptions(solver);
280 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
281 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
290 Tag th_error_ind, th_order;
294 std::vector<Range> refinement_levels;
297 double err_indic = 0;
299 int order, new_order;
301 new_order =
order + 1;
304 refined_ents.insert(ent);
307 moab::Interface::UNION);
308 refined_ents.merge(adj);
309 refinement_levels[new_order -
baseOrder].merge(refined_ents);
314 for (
int ll = 1; ll < refinement_levels.size(); ll++) {
316 refinement_levels[ll]);
362 auto det_ptr = boost::make_shared<VectorDouble>();
363 auto jac_ptr = boost::make_shared<MatrixDouble>();
364 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
408 <<
"Global error L2 norm: " << std::sqrt(array[0]);
410 <<
"Global error H1 seminorm: " << std::sqrt(array[1]);
412 <<
"Global error indicator: " << std::sqrt(array[2]);
414 <<
"Total number of elements: " << (int)array[3];
420 std::vector<Tag> tag_handles;
421 tag_handles.resize(4);
429 ParallelComm *pcomm =
434 tag_handles.push_back(pcomm->part_tag());
435 std::ostringstream strm;
436 strm <<
"error_" << iter_num <<
".h5m";
438 "PARALLEL=WRITE_PART", 0, 0,
439 tag_handles.data(), tag_handles.size());
453 boost::make_shared<PostProcEle>(
mField);
455 auto det_ptr = boost::make_shared<VectorDouble>();
456 auto jac_ptr = boost::make_shared<MatrixDouble>();
457 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
459 post_proc_fe->getOpPtrVector().push_back(
461 post_proc_fe->getOpPtrVector().push_back(
464 post_proc_fe->getOpPtrVector().push_back(
466 post_proc_fe->getOpPtrVector().push_back(
469 auto u_ptr = boost::make_shared<VectorDouble>();
470 auto flux_ptr = boost::make_shared<MatrixDouble>();
471 post_proc_fe->getOpPtrVector().push_back(
473 post_proc_fe->getOpPtrVector().push_back(
478 post_proc_fe->getOpPtrVector().push_back(
480 new OpPPMap(post_proc_fe->getPostProcMesh(),
481 post_proc_fe->getMapGaussPts(),
495 pipeline_mng->getDomainRhsFE() = post_proc_fe;
496 CHKERR pipeline_mng->loopFiniteElements();
498 std::ostringstream strm;
499 strm <<
"out_" << iter_num <<
".h5m";
500 CHKERR post_proc_fe->writeFile(strm.str().c_str());
509 const int nb_integration_pts =
getGaussPts().size2();
513 auto t_val_grad = getFTensor1FromMat<2>(*(
commonDataPtr->approxValsGrad));
514 auto t_flux = getFTensor1FromMat<3>(*(
commonDataPtr->approxFlux));
521 double error_ind = 0;
522 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
523 const double alpha = t_w * area;
526 t_coords(0), t_coords(1), t_coords(2));
527 error_l2 += alpha *
sqr(diff);
530 t_coords(0), t_coords(1), t_coords(2));
531 auto t_fun_grad = getFTensor1FromArray<2, 2>(vec);
532 t_diff(
i) = t_val_grad(
i) - t_fun_grad(
i);
533 error_h1 += alpha * t_diff(
i) * t_diff(
i);
535 t_diff(
i) = t_val_grad(
i) - t_flux(
i);
536 error_ind += alpha * t_diff(
i) * t_diff(
i);
546 Tag th_error_l2, th_error_h1, th_error_ind;
558 constexpr std::array<int, 4> indices = {
561 std::array<double, 4> values;
562 values[0] = error_l2;
563 values[1] = error_h1;
564 values[2] = error_ind;
572int main(
int argc,
char *argv[]) {
578 auto core_log = logging::core::get();
586 DMType dm_name =
"DMMOFEM";
591 moab::Core mb_instance;
592 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_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 >::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.
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
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 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
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]
double errorIndicatorIntegral
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name, DataType type, Tag &tag_handle)
[Source function]
MoFEMErrorCode readMesh()
[Run programme]
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.
Make Hdiv space from Hcurl space in 2d.
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
Modify integration weights on face to take in account higher-order geometry.
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]