Testing projection child and parent.
static char help[] =
"...\n\n";
};
};
double operator()(const double x, const double y, const double z) {
return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) +
pow(y, 4);
}
};
private:
};
template <
int FIELD_DIM>
struct OpError;
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
if (
const size_t nb_dofs = data.
getIndices().size()) {
nf.clear();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
t_coords(2));
error +=
alpha * pow(diff, 2);
for (
size_t r = 0;
r != nb_dofs; ++
r) {
nf[
r] +=
alpha * t_row_base * diff;
++t_row_base;
}
++t_w;
++t_val;
++t_coords;
}
}
}
};
template <typename ELE_OP, typename PARENT_ELE>
PARENT_ELE parent_fe(this->getPtrFE()->mField);
<< "parent_coords in op "
<<
static_cast<ELE_OP *
>(op_ptr)->getCoordsAtGaussPts();
parent_coords =
static_cast<ELE_OP *
>(op_ptr)->getCoordsAtGaussPts();
};
parent_fe.getOpPtrVector().push_back(op);
MOFEM_LOG(
"SELF", Sev::noisy) <<
"fe name " << this->getFEName();
CHKERR this->loopParent(this->getFEName(), &parent_fe);
MOFEM_LOG(
"SELF", Sev::noisy) <<
"parent_coords " << parent_coords;
MOFEM_LOG(
"SELF", Sev::noisy) <<
"child_coords " << child_coords;
child_coords -= parent_coords;
MOFEM_LOG(
"SELF", Sev::noisy) <<
"Corrds diffs" << child_coords;
for (auto d : child_coords.data())
"Parent and child global coords at integration points are "
"diffrent norm = %3.2e",
}
};
const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
};
}
auto refine_mesh = [&](auto bit_level1) {
CHKERR moab.get_entities_by_type(*meshset_level0_ptr, MBEDGE,
edges_to_refine);
int ii = 0;
for (Range::iterator eit = edges_to_refine.begin();
eit != edges_to_refine.end(); eit++, ii++) {
int numb = ii % 2;
if (numb == 0) {
CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*eit, 1);
}
}
CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
CHKERR refine->refineTets(*meshset_level0_ptr, bit_level1,
VERBOSE);
CHKERR refine->refineTris(*meshset_level0_ptr, bit_level1,
VERBOSE);
} else {
"Dimension not handled by test");
}
};
bit_level1.set(1);
CHKERR refine_mesh(bit_level1);
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto test_bit_parent = [](
FEMethod *fe_ptr) {
const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
};
parent_op_lhs->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
MOFEM_LOG(
"SELF", Sev::noisy) <<
"LHS Pipeline FE";
domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
CHKERR domain_op->loopChildren(domain_op->getFEName(),
} else {
}
};
parent_op_rhs->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
MOFEM_LOG(
"SELF", Sev::noisy) <<
"RHS Pipeline FE";
domain_op->getFEMethod()->numeredEntFiniteElementPtr->getBitRefLevel();
CHKERR domain_op->loopChildren(domain_op->getFEName(),
}
};
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
auto refine_mesh = [&]() {
*meshset_level1_ptr);
bit_level1,
BitRefLevel().set(), MBEDGE, edges_to_refine);
CHKERR refine->addVerticesInTheMiddleOfEdges(edges_to_refine, bit_level2,
CHKERR refine->refineTets(*meshset_level1_ptr, bit_level2,
VERBOSE);
CHKERR refine->refineTris(*meshset_level1_ptr, bit_level2,
VERBOSE);
} else {
"Dimension not handled by test");
}
CHKERR moab.get_entities_by_type(0, MBENTITYSET, meshsets,
true);
for (
auto m : meshsets) {
->updateMeshsetByEntitiesChildren(
m, bit_level2,
m, MBMAXTYPE,
false);
}
};
bit_level0 | bit_level1);
auto project_data = [&]() {
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto test_bit_ref = [](
FEMethod *fe_ptr) {
const auto &
bit = fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
};
auto field_vals_ptr = boost::make_shared<VectorDouble>();
auto domainParentRhs = boost::make_shared<DomainParentEle>(
mField);
domainParentRhs->getOpPtrVector().push_back(
new OpRunParent(domainParentRhs, bit_level2, bit_level2,
};
}
boost::function<
bool(
FEMethod *fe_method_ptr)> test_bit) {
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
auto common_data_ptr = boost::make_shared<CommonData>();
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
common_data_ptr->approxVals));
CHKERR VecZeroEntries(common_data_ptr->L2Vec);
CHKERR VecZeroEntries(common_data_ptr->resVec);
CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
CHKERR VecAssemblyBegin(common_data_ptr->resVec);
CHKERR VecAssemblyEnd(common_data_ptr->resVec);
double nrm2;
CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
const double *array;
CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
MOFEM_LOG_C(
"WORLD", Sev::inform,
"Error %6.4e Vec norm %6.4e\n",
std::sqrt(array[0]), nrm2);
CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
constexpr double eps = 1e-8;
"Not converged solution err = %6.4e", nrm2);
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
#define MOFEM_LOG_C(channel, severity, format,...)
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
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]
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#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.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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.
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 > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< 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
constexpr char FIELD_NAME[]
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
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.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
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.
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
static ApproxFieldFunction< FIELD_DIM > approxFunction
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.
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.
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.
intrusive_ptr for managing petsc objects
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)