12static char help[] =
"...\n\n";
39 auto operator()(
const double x,
const double y,
const double z) {
40 return x * x + y * y + x * y * y + x * x * y;
49 auto operator()(
const double x,
const double y,
const double z) {
53 2 * y + 2 * x * y + x * x};
91template <
typename PARENT_FE>
93 boost::shared_ptr<FEMethod> &fe_top,
94 ForcesAndSourcesCore::UserDataOperator::OpType op,
104 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
106 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
112 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
113 new PARENT_FE(m_field));
118 fe_ptr_current->getOpPtrVector(), {H1});
123 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
131 parent_fe_pt->getOpPtrVector().push_back(
135 H1, op, fe_ptr_current,
146 parent_fe_pt->getOpPtrVector().push_back(
161 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
174 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
215 template <
int FIELD_DIM>
struct OpError;
226 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
231 if (
const size_t nb_dofs = data.
getIndices().size()) {
235 const int nb_integration_pts =
getGaussPts().size2();
239 getFTensor1FromMat<SPACE_DIM>(*(
commonDataPtr->divApproxVals));
249 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
251 const double alpha = t_w * volume;
257 t_grad_diff(
i) -= t_grad_val(
i);
259 error += alpha * (pow(diff, 2) + t_grad_diff(
i) * t_grad_diff(
i));
261 for (
size_t r = 0; r != nb_dofs; ++r) {
262 nf[r] += alpha * t_row_base * diff;
290 const int nb_integration_pts = getGaussPts().size2();
291 auto t_w = getFTensor0IntegrationWeight();
293 auto t_coords = getFTensor1CoordsAtGaussPts();
295 const double volume = getMeasure();
298 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
300 const double alpha = t_w * volume;
303 error2 += alpha * (pow(diff, 2));
310 MOFEM_LOG(
"SELF", Sev::verbose) <<
"Boundary error " << sqrt(error2);
312 constexpr double eps = 1e-8;
313 if (sqrt(error2) >
eps)
315 "Error on boundary = %6.4e", sqrt(error2));
337 ParallelComm *pcomm =
354 CHKERR skin.find_skin(0, level0_ents,
false, level0_skin);
355 CHKERR pcomm->filter_pstatus(level0_skin,
356 PSTATUS_SHARED | PSTATUS_MULTISHARED,
357 PSTATUS_NOT, -1,
nullptr);
359 auto refine_mesh = [&](
auto l) {
372 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr,
SPACE_DIM, els);
381 ele_to_refine.insert(
t);
382 std::vector<EntityHandle> adj_edges;
385 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
392 CHKERR skin.find_skin(0, els,
false, level_skin);
393 CHKERR pcomm->filter_pstatus(level_skin,
394 PSTATUS_SHARED | PSTATUS_MULTISHARED,
395 PSTATUS_NOT, -1,
nullptr);
396 level_skin = subtract(level_skin, level0_skin);
399 adj, moab::Interface::UNION);
400 els = subtract(els, adj);
401 ele_to_refine.merge(els);
404 els,
SPACE_DIM - 1,
false, adj_edges, moab::Interface::UNION);
405 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
408 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
bit(
l),
419 (boost::lexical_cast<std::string>(
l) +
"_ref_mesh.vtk").c_str(),
"VTK",
423 (boost::lexical_cast<std::string>(
l) +
"_only_ref_mesh.vtk").c_str(),
429 auto mark_skins = [&](
auto l,
auto m) {
435 CHKERR skin.find_skin(0, ents,
false, level_skin);
436 CHKERR pcomm->filter_pstatus(level_skin,
437 PSTATUS_SHARED | PSTATUS_MULTISHARED,
438 PSTATUS_NOT, -1,
nullptr);
439 level_skin = subtract(level_skin, level0_skin);
441 moab::Interface::UNION);
471 constexpr int order = 3;
501 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
521 field_op_row->doWorkRhsHook = [](
DataOperator *op_ptr,
int side,
525 if (type == MBENTITYSET) {
528 <<
"ROW: side/type: " << side <<
"/" << CN::EntityTypeName(type)
531 << data.getIndices() <<
" nb base functions " << data.getN().size2()
532 <<
" nb base functions integration points " << data.getN().size1();
534 for (
auto &field_ent : data.getFieldEntities()) {
536 <<
"\t" << CN::EntityTypeName(field_ent->getEntType());
560 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
563 CHKERR KSPSetFromOptions(solver);
571 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
572 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
585 auto rule = [](int, int,
int p) ->
int {
return 2 *
p + 1; };
591 auto common_data_ptr = boost::make_shared<CommonData>();
595 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
596 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
610 common_data_ptr->approxVals));
621 common_data_ptr->approxVals));
625 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
626 CHKERR VecZeroEntries(common_data_ptr->resVec);
630 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
631 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
632 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
633 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
635 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
637 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
638 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Error %6.4e Vec norm %6.4e\n",
639 std::sqrt(array[0]), nrm2);
641 constexpr double eps = 1e-8;
644 "Not converged solution err = %6.4e", nrm2);
645 if (std::sqrt(array[0]) >
eps)
647 "Error in approximation err = %6.4e", std::sqrt(array[0]));
649 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
661 auto rule = [](int, int,
int p) ->
int {
return -1; };
684 auto approx_vals = boost::make_shared<VectorDouble>();
688 double def_val[] = {0};
689 CHKERR moab.tag_get_handle(
"FIELD", 1, MB_TYPE_DOUBLE,
th,
690 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
692 field_op_row->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
696 if (type == MBVERTEX) {
701 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
702 if (nb_gauss_pts != 3)
704 "Should be three guass pts.");
705 auto conn = op_ptr->getConn();
706 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
707 const double v = t_field;
708 CHKERR moab.tag_set_data(
th, &conn[gg], 1, &
v);
730int main(
int argc,
char *argv[]) {
738 DMType dm_name =
"DMMOFEM";
743 moab::Core mb_instance;
744 moab::Interface &moab = mb_instance;
#define MOFEM_LOG_C(channel, severity, format,...)
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
#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 ...
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.
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 addBitRefLevel(const Range &ents, const BitRefLevel &bit, int verb=QUIET) const
add bit ref level to ref entity
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 updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
auto marker
set bit to marker
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.
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
constexpr char FIELD_NAME[]
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.
const 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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI 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 temporary meshset.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr double t
plate stiffness
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()
Add operators pushing bases from local to physical configuration.
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.
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
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 to child.
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.
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.