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,
97 auto jac_ptr = boost::make_shared<MatrixDouble>();
98 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
99 auto det_ptr = boost::make_shared<VectorDouble>();
108 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
110 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
116 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
117 new PARENT_FE(m_field));
122 fe_ptr_current->getOpPtrVector(), {H1});
127 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
135 parent_fe_pt->getOpPtrVector().push_back(
139 H1, op, fe_ptr_current,
150 parent_fe_pt->getOpPtrVector().push_back(
165 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
178 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
219 template <
int FIELD_DIM>
struct OpError;
230 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
235 if (
const size_t nb_dofs = data.
getIndices().size()) {
239 const int nb_integration_pts =
getGaussPts().size2();
243 getFTensor1FromMat<SPACE_DIM>(*(
commonDataPtr->divApproxVals));
253 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
255 const double alpha = t_w * volume;
261 t_grad_diff(
i) -= t_grad_val(
i);
263 error += alpha * (pow(diff, 2) + t_grad_diff(
i) * t_grad_diff(
i));
265 for (
size_t r = 0; r != nb_dofs; ++r) {
266 nf[r] += alpha * t_row_base * diff;
294 const int nb_integration_pts = getGaussPts().size2();
295 auto t_w = getFTensor0IntegrationWeight();
297 auto t_coords = getFTensor1CoordsAtGaussPts();
299 const double volume = getMeasure();
302 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
304 const double alpha = t_w * volume;
307 error2 += alpha * (pow(diff, 2));
314 MOFEM_LOG(
"SELF", Sev::verbose) <<
"Boundary error " << sqrt(error2);
316 constexpr double eps = 1e-8;
317 if (sqrt(error2) >
eps)
319 "Error on boundary = %6.4e", sqrt(error2));
341 ParallelComm *pcomm =
358 CHKERR skin.find_skin(0, level0_ents,
false, level0_skin);
359 CHKERR pcomm->filter_pstatus(level0_skin,
360 PSTATUS_SHARED | PSTATUS_MULTISHARED,
361 PSTATUS_NOT, -1,
nullptr);
363 auto refine_mesh = [&](
auto l) {
376 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr,
SPACE_DIM, els);
385 ele_to_refine.insert(
t);
386 std::vector<EntityHandle> adj_edges;
389 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
396 CHKERR skin.find_skin(0, els,
false, level_skin);
397 CHKERR pcomm->filter_pstatus(level_skin,
398 PSTATUS_SHARED | PSTATUS_MULTISHARED,
399 PSTATUS_NOT, -1,
nullptr);
400 level_skin = subtract(level_skin, level0_skin);
403 adj, moab::Interface::UNION);
404 els = subtract(els, adj);
405 ele_to_refine.merge(els);
408 els,
SPACE_DIM - 1,
false, adj_edges, moab::Interface::UNION);
409 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
412 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
bit(
l),
423 (boost::lexical_cast<std::string>(
l) +
"_ref_mesh.vtk").c_str(),
"VTK",
427 (boost::lexical_cast<std::string>(
l) +
"_only_ref_mesh.vtk").c_str(),
433 auto mark_skins = [&](
auto l,
auto m) {
439 CHKERR skin.find_skin(0, ents,
false, level_skin);
440 CHKERR pcomm->filter_pstatus(level_skin,
441 PSTATUS_SHARED | PSTATUS_MULTISHARED,
442 PSTATUS_NOT, -1,
nullptr);
443 level_skin = subtract(level_skin, level0_skin);
445 moab::Interface::UNION);
475 constexpr int order = 3;
506 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
526 field_op_row->doWorkRhsHook = [](
DataOperator *op_ptr,
int side,
530 if (type == MBENTITYSET) {
533 <<
"ROW: side/type: " << side <<
"/" << CN::EntityTypeName(type)
536 << data.getIndices() <<
" nb base functions " << data.getN().size2()
537 <<
" nb base functions integration points " << data.getN().size1();
539 for (
auto &field_ent : data.getFieldEntities()) {
541 <<
"\t" << CN::EntityTypeName(field_ent->getEntType());
565 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
568 CHKERR KSPSetFromOptions(solver);
576 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
577 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
590 auto rule = [](int, int,
int p) ->
int {
return 2 *
p + 1; };
596 auto common_data_ptr = boost::make_shared<CommonData>();
600 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
601 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
615 common_data_ptr->approxVals));
626 common_data_ptr->approxVals));
630 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
631 CHKERR VecZeroEntries(common_data_ptr->resVec);
635 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
636 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
637 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
638 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
640 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
642 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
643 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Error %6.4e Vec norm %6.4e\n",
644 std::sqrt(array[0]), nrm2);
646 constexpr double eps = 1e-8;
649 "Not converged solution err = %6.4e", nrm2);
650 if (std::sqrt(array[0]) >
eps)
652 "Error in approximation err = %6.4e", std::sqrt(array[0]));
654 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
666 auto rule = [](int, int,
int p) ->
int {
return -1; };
689 auto approx_vals = boost::make_shared<VectorDouble>();
693 double def_val[] = {0};
694 CHKERR moab.tag_get_handle(
"FIELD", 1, MB_TYPE_DOUBLE,
th,
695 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
697 field_op_row->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
701 if (type == MBVERTEX) {
706 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
707 if (nb_gauss_pts != 3)
709 "Should be three guass pts.");
710 auto conn = op_ptr->getConn();
711 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
712 const double v = t_field;
713 CHKERR moab.tag_set_data(
th, &conn[gg], 1, &
v);
735int main(
int argc,
char *argv[]) {
743 DMType dm_name =
"DMMOFEM";
748 moab::Core mb_instance;
749 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 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 > & 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.
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
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 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.