static char help[] =
"...\n\n";
GAUSS>::OpMixDivTimesScalar<2>;
inline double sqr(
double x) {
return x * x; }
private:
const double z) {
return exp(-100. * (
sqr(x) +
sqr(y))) * cos(M_PI * x) * cos(M_PI * 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;
switch (type) {
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);
}
}
};
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr,
}
};
};
}
}
1);
PETSC_NULL);
PETSC_NULL);
false);
Tag th_order;
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
}
PetscInt ghosts[4] = {0, 1, 2, 3};
else
commonDataPtr->approxValsGrad = boost::make_shared<MatrixDouble>();
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto unity = []() { return 1; };
new OpHdivU(
"FLUX",
"U", unity,
true));
auto source = [&](
const double x,
const double y,
const double z) {
};
}
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
Tag th_error_ind, th_order;
std::vector<Range> refinement_levels;
double err_indic = 0;
refined_ents.insert(ent);
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]);
}
}
}
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
SCATTER_FORWARD);
SCATTER_FORWARD);
SCATTER_REVERSE);
SCATTER_REVERSE);
SCATTER_FORWARD);
SCATTER_FORWARD);
const double *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];
std::vector<Tag> tag_handles;
tag_handles.resize(4);
tag_handles[1]);
tag_handles[2]);
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());
}
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(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<VectorDouble>();
auto flux_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
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{}
)
);
std::ostringstream strm;
strm << "out_" << iter_num << ".h5m";
CHKERR post_proc_fe->writeFile(strm.str().c_str());
}
auto t_val_grad = getFTensor1FromMat<2>(*(
commonDataPtr->approxValsGrad));
auto t_flux = getFTensor1FromMat<3>(*(
commonDataPtr->approxFlux));
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);
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.;
ADD_VALUES);
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
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.
PetscErrorCode DMSetUp_MoFEM(DM dm)
auto smartCreateDMVector(DM dm)
Get smart vector from 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.
#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.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
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]
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]
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
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.
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]