26static char help[] =
"...\n\n";
60 double operator()(
const double x,
const double y,
const double z) {
61 return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) +
96 template <
int FIELD_DIM>
struct OpError;
104 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
109 if (
const size_t nb_dofs = data.
getIndices().size()) {
111 const int nb_integration_pts =
getGaussPts().size2();
124 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
126 const double alpha = t_w * volume;
129 error +=
alpha * pow(diff, 2);
131 for (
size_t r = 0;
r != nb_dofs; ++
r) {
132 nf[
r] +=
alpha * t_row_base * diff;
150template <
typename ELE_OP,
typename PARENT_ELE>
160 PARENT_ELE parent_fe(this->getPtrFE()->mField);
168 <<
"parent_coords in op "
169 <<
static_cast<ELE_OP *
>(op_ptr)->getCoordsAtGaussPts();
171 parent_coords =
static_cast<ELE_OP *
>(op_ptr)->getCoordsAtGaussPts();
174 parent_fe.getOpPtrVector().push_back(op);
176 MOFEM_LOG(
"SELF", Sev::noisy) <<
"fe name " << this->getFEName();
177 CHKERR this->loopParent(this->getFEName(), &parent_fe);
178 MOFEM_LOG(
"SELF", Sev::noisy) <<
"parent_coords " << parent_coords;
180 MatrixDouble child_coords = this->getCoordsAtGaussPts();
181 MOFEM_LOG(
"SELF", Sev::noisy) <<
"child_coords " << child_coords;
183 child_coords -= parent_coords;
185 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Corrds diffs" << child_coords;
188 for (
auto d : child_coords.data())
193 "Parent and child global coords at integration points are "
194 "diffrent norm = %3.2e",
210 const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
235 auto refine_mesh = [&](
auto bit_level1) {
242 bit_level0,
BitRefLevel().set(), *meshset_level0_ptr);
246 Range edges_to_refine;
247 CHKERR moab.get_entities_by_type(*meshset_level0_ptr, MBEDGE,
250 for (Range::iterator eit = edges_to_refine.begin();
251 eit != edges_to_refine.end(); eit++, ii++) {
254 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*eit, 1);
257 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
260 CHKERR refine->refineTets(*meshset_level0_ptr, bit_level1,
VERBOSE);
262 CHKERR refine->refineTris(*meshset_level0_ptr, bit_level1,
VERBOSE);
265 "Dimension not handled by test");
273 CHKERR refine_mesh(bit_level1);
289 constexpr int order = 4;
308 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
313 auto test_bit_parent = [](
FEMethod *fe_ptr) {
314 const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
336 parent_op_lhs->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
339 auto domain_op =
static_cast<DomainEleOp *
>(op_ptr);
342 MOFEM_LOG(
"SELF", Sev::noisy) <<
"LHS Pipeline FE";
348 domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
350 CHKERR domain_op->loopChildren(domain_op->getFEName(),
360 parent_op_rhs->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
363 auto domain_op =
static_cast<DomainEleOp *
>(op_ptr);
366 MOFEM_LOG(
"SELF", Sev::noisy) <<
"RHS Pipeline FE";
372 domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
374 CHKERR domain_op->loopChildren(domain_op->getFEName(),
395 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
398 CHKERR KSPSetFromOptions(solver);
406 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
407 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
422 auto refine_mesh = [&]() {
430 *meshset_level1_ptr);
434 Range edges_to_refine;
436 bit_level1,
BitRefLevel().set(), MBEDGE, edges_to_refine);
438 CHKERR refine->addVerticesInTheMiddleOfEdges(edges_to_refine, bit_level2,
441 CHKERR refine->refineTets(*meshset_level1_ptr, bit_level2,
VERBOSE);
443 CHKERR refine->refineTris(*meshset_level1_ptr, bit_level2,
VERBOSE);
446 "Dimension not handled by test");
450 CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets,
true);
451 for (
auto m : meshsets) {
453 ->updateMeshsetByEntitiesChildren(
m, bit_level2,
m, MBMAXTYPE,
false);
468 bit_level0 | bit_level1);
470 auto project_data = [&]() {
479 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
485 auto test_bit_ref = [](
FEMethod *fe_ptr) {
486 const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
487 MOFEM_LOG(
"SELF", Sev::noisy) <<
"ref : " <<
bit <<
" " <<
bit.test(2);
496 auto field_vals_ptr = boost::make_shared<VectorDouble>();
498 auto domainParentRhs = boost::make_shared<DomainParentEle>(
mField);
499 domainParentRhs->getOpPtrVector().push_back(
505 new OpRunParent(domainParentRhs, bit_level2, bit_level2,
506 domainParentRhs, bit_level2,
BitRefLevel().set()));
538 boost::function<
bool(
FEMethod *fe_method_ptr)> test_bit) {
545 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
549 auto common_data_ptr = boost::make_shared<CommonData>();
553 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
557 common_data_ptr->approxVals));
561 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
562 CHKERR VecZeroEntries(common_data_ptr->resVec);
566 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
567 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
568 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
569 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
571 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
573 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
574 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Error %6.4e Vec norm %6.4e\n",
575 std::sqrt(array[0]), nrm2);
576 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
578 constexpr double eps = 1e-8;
581 "Not converged solution err = %6.4e", nrm2);
586int main(
int argc,
char *argv[]) {
594 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
boost::shared_ptr< DomainEle > domainChildRhs
boost::shared_ptr< DomainEle > domainChildLhs
[Set up problem]
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
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 MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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.
#define MOFEM_LOG(channel, severity)
Log.
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
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
FTensor::Index< 'm', 3 > m
double 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
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
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.
Operator to evaluate errors.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
static ApproxFieldFunction< FIELD_DIM > approxFunction
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode refineResults()
[Solve]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Run programme]
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
@ 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
Mesh refinement interface.
Get value at integration points for scalar field.
Operator to execute finite element instance on parent element. This operator is typically used to pro...
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 reSetUp(bool only_dm=false)
Rebuild internal MoFEM data structures.
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 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.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getProblemName() const
Get the Problem Name.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)