26static char help[] =
"...\n\n";
66 auto operator()(
const double x,
const double y,
const double z) {
67 return x * x + y * y + x * y * y + x * x * y;
76 auto operator()(
const double x,
const double y,
const double z) {
80 2 * y + 2 * x * y + x * x};
118template <
typename PARENT_FE>
120 boost::shared_ptr<FEMethod> &fe_top,
121 ForcesAndSourcesCore::UserDataOperator::OpType op,
124 auto jac_ptr = boost::make_shared<MatrixDouble>();
125 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
126 auto det_ptr = boost::make_shared<VectorDouble>();
132 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
134 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
137 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
138 new PARENT_FE(m_field));
141 fe_ptr_current->getOpPtrVector().push_back(
143 fe_ptr_current->getOpPtrVector().push_back(
145 fe_ptr_current->getOpPtrVector().push_back(
151 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
157 parent_fe_pt->getOpPtrVector().push_back(
161 H1, op, fe_ptr_current,
171 parent_fe_pt->getOpPtrVector().push_back(
186 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
199 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
240 template <
int FIELD_DIM>
struct OpError;
251 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
256 if (
const size_t nb_dofs = data.
getIndices().size()) {
260 const int nb_integration_pts =
getGaussPts().size2();
264 getFTensor1FromMat<SPACE_DIM>(*(
commonDataPtr->divApproxVals));
274 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
276 const double alpha = t_w * volume;
282 t_grad_diff(
i) -= t_grad_val(
i);
284 error +=
alpha * (pow(diff, 2) + t_grad_diff(
i) * t_grad_diff(
i));
286 for (
size_t r = 0;
r != nb_dofs; ++
r) {
287 nf[
r] +=
alpha * t_row_base * diff;
315 const int nb_integration_pts = getGaussPts().size2();
316 auto t_w = getFTensor0IntegrationWeight();
318 auto t_coords = getFTensor1CoordsAtGaussPts();
320 const double volume = getMeasure();
323 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
325 const double alpha = t_w * volume;
328 error2 +=
alpha * (pow(diff, 2));
335 MOFEM_LOG(
"SELF", Sev::verbose) <<
"Boundary error " << sqrt(error2);
337 constexpr double eps = 1e-8;
338 if (sqrt(error2) >
eps)
340 "Error on boundary = %6.4e", sqrt(error2));
363 ParallelComm *pcomm =
380 CHKERR skin.find_skin(0, level0_ents,
false, level0_skin);
381 CHKERR pcomm->filter_pstatus(level0_skin,
382 PSTATUS_SHARED | PSTATUS_MULTISHARED,
383 PSTATUS_NOT, -1,
nullptr);
385 auto refine_mesh = [&](
auto l) {
398 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr,
SPACE_DIM, els);
407 ele_to_refine.insert(
t);
408 std::vector<EntityHandle> adj_edges;
411 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
418 CHKERR skin.find_skin(0, els,
false, level_skin);
419 CHKERR pcomm->filter_pstatus(level_skin,
420 PSTATUS_SHARED | PSTATUS_MULTISHARED,
421 PSTATUS_NOT, -1,
nullptr);
422 level_skin = subtract(level_skin, level0_skin);
425 adj, moab::Interface::UNION);
426 els = subtract(els, adj);
427 ele_to_refine.merge(els);
430 els,
SPACE_DIM - 1,
false, adj_edges, moab::Interface::UNION);
431 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
434 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
bit(
l),
445 (boost::lexical_cast<std::string>(
l) +
"_ref_mesh.vtk").c_str(),
"VTK",
449 (boost::lexical_cast<std::string>(
l) +
"_only_ref_mesh.vtk").c_str(),
455 auto mark_skins = [&](
auto l,
auto m) {
461 CHKERR skin.find_skin(0, ents,
false, level_skin);
462 CHKERR pcomm->filter_pstatus(level_skin,
463 PSTATUS_SHARED | PSTATUS_MULTISHARED,
464 PSTATUS_NOT, -1,
nullptr);
465 level_skin = subtract(level_skin, level0_skin);
467 moab::Interface::UNION);
497 constexpr int order = 3;
532 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
552 field_op_row->doWorkRhsHook = [](
DataOperator *op_ptr,
int side,
556 if (
type == MBENTITYSET) {
559 <<
"ROW: side/type: " << side <<
"/" << CN::EntityTypeName(
type)
562 << data.getIndices() <<
" nb base functions " << data.getN().size2()
563 <<
" nb base functions integration points " << data.getN().size1();
565 auto get_bool = [](
auto fe,
auto bit) {
566 return (
bit & fe->getBitRefLevel()).any();
569 for (
auto &field_ent : data.getFieldEntities()) {
571 <<
"\t" << CN::EntityTypeName(field_ent->getEntType());
577 set_parent_dofs<DomainParentEle>(
596 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
599 CHKERR KSPSetFromOptions(solver);
607 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
608 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
621 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
627 auto common_data_ptr = boost::make_shared<CommonData>();
631 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
632 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
634 auto jac_ptr = boost::make_shared<MatrixDouble>();
635 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
636 auto det_ptr = boost::make_shared<VectorDouble>();
655 common_data_ptr->approxVals));
666 common_data_ptr->approxVals));
670 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
671 CHKERR VecZeroEntries(common_data_ptr->resVec);
675 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
676 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
677 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
678 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
680 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
682 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
683 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Error %6.4e Vec norm %6.4e\n",
684 std::sqrt(array[0]), nrm2);
686 constexpr double eps = 1e-8;
689 "Not converged solution err = %6.4e", nrm2);
690 if (std::sqrt(array[0]) >
eps)
692 "Error in approximation err = %6.4e", std::sqrt(array[0]));
694 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
706 auto rule = [](
int,
int,
int p) ->
int {
return -1; };
729 auto approx_vals = boost::make_shared<VectorDouble>();
733 double def_val[] = {0};
734 CHKERR moab.tag_get_handle(
"FIELD", 1, MB_TYPE_DOUBLE,
th,
735 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
737 field_op_row->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
741 if (
type == MBVERTEX) {
746 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
747 if (nb_gauss_pts != 3)
749 "Should be three guass pts.");
750 auto conn = op_ptr->getConn();
751 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
752 const double v = t_field;
753 CHKERR moab.tag_set_data(
th, &conn[gg], 1, &
v);
775int main(
int argc,
char *argv[]) {
783 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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 smartCreateDMVector(DM dm)
Get smart vector from DM.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
MoFEMErrorCode updateMeshsetByEntitiesChildren(const EntityHandle parent, const BitRefLevel &parent_bit, const BitRefLevel &parent_mask, const BitRefLevel &child_bit, const BitRefLevel &child_mask, const EntityHandle child, EntityType child_type, const bool recursive=false, int verb=0)
Get child entities form meshset containing parent entities.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel bit, int verb=QUIET) const
add bit ref level to ref entity
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childresn.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
auto marker
set bit to marker
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
constexpr int nb_ref_levels
Three levels of refinement.
constexpr char FIELD_NAME[]
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
evaluate source, i.e. rhs vector
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
evaluate mass matrix
auto set_parent_dofs(MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > &fe_top, ForcesAndSourcesCore::UserDataOperator::OpType op, int verbosity, LogManager::SeverityLevel sev)
set levels of projection operators, which project field data from parent entities,...
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< 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 get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
DeprecatedCoreInterface Interface
constexpr auto VecSetValues
const double D
diffusivity
const double r
rate factor
constexpr double t
plate stiffness
FTensor::Index< 'm', 3 > m
auto operator()(const double x, const double y, const double z)
auto operator()(const double x, const double y, const double z)
Collected data use d by operator to evaluate errors for the test.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< MatrixDouble > divApproxVals
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator to evaluate errors.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpErrorSkel(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
static ApproxFieldFunction< FIELD_DIM > approxFunction
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
add field, and set up problem
MoFEMErrorCode readMesh()
red mesh and randomly refine three times
MoFEMErrorCode assembleSystem()
static ApproxFieldFunctionDerivative< FIELD_DIM > divApproxFunction
MoFEM::Interface & mField
MoFEMErrorCode printResults()
[Check results]
MoFEMErrorCode runProblem()
MoFEMErrorCode writeBitLevelByDim(const BitRefLevel bit, const BitRefLevel mask, const int dim, const char *file_name, const char *file_type, const char *options, const bool check_for_empty=true) const
Write bit ref level to file.
MoFEMErrorCode filterEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, Range &ents, int verb=QUIET) const
filter entities by bit ref level
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.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
default operator for EDGE element
Base face element used to integrate on skeleton.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
default operator for TRI element
Base face element used to integrate on skeleton.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
Mesh refinement interface.
Operator to project base functions from parent entity.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
int getDim() const
Get the problem dimension.
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.
auto & getBitAdjEnt()
bit ref level for parent
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 addBoundaryField(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 boundary.
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevel.
EntityHandle & getBoundaryMeshSet()
Get the BoundaryMeshSet object.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
bool & getParentAdjacencies()
Get the addParentAdjacencies.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.