Testing approximation with hanging nodes.
static char help[] =
"...\n\n";
auto operator()(const double x, const double y, const double z) {
return x * x + y * y + x * y * y + x * x * y;
}
};
auto operator()(const double x, const double y, const double z) {
2 * y + 2 * x * y + x * x};
}
};
};
template <typename PARENT_FE>
boost::shared_ptr<FEMethod> &fe_top,
ForcesAndSourcesCore::UserDataOperator::OpType op,
boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
add_parent_level =
[&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
if (level > 0) {
auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
new PARENT_FE(m_field));
fe_ptr_current->getOpPtrVector(), {H1});
}
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
parent_fe_pt->getOpPtrVector().push_back(
verbosity, sev));
} else {
parent_fe_pt->getOpPtrVector().push_back(
verbosity, sev));
}
}
};
add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
};
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
private:
};
template <
int FIELD_DIM>
struct OpError;
template <int FIELD_DIM> struct OpErrorSkel;
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
if (
const size_t nb_dofs = data.
getIndices().size()) {
auto t_grad_val =
nf.clear();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
t_coords(2));
auto t_grad_diff =
t_grad_diff(
i) -= t_grad_val(
i);
error +=
alpha * (pow(diff, 2) + t_grad_diff(
i) * t_grad_diff(
i));
for (
size_t r = 0;
r != nb_dofs; ++
r) {
nf[
r] +=
alpha * t_row_base * diff;
++t_row_base;
}
++t_w;
++t_val;
++t_grad_val;
++t_coords;
}
}
}
};
boost::shared_ptr<CommonData> commonDataPtr;
OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
auto t_coords = getFTensor1CoordsAtGaussPts();
const double volume = getMeasure();
double error2 = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
t_coords(2));
error2 +=
alpha * (pow(diff, 2));
++t_w;
++t_val;
++t_coords;
}
MOFEM_LOG(
"SELF", Sev::verbose) <<
"Boundary error " << sqrt(error2);
constexpr double eps = 1e-8;
"Error on boundary = %6.4e", sqrt(error2));
}
};
}
ParallelComm *pcomm =
CHKERR skin.find_skin(0, level0_ents,
false, level0_skin);
CHKERR pcomm->filter_pstatus(level0_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
auto refine_mesh = [&](
auto l) {
CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr,
SPACE_DIM, els);
int ii = 0;
if ((ii % 2)) {
std::vector<EntityHandle> adj_edges;
adj_edges);
CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
adj_edges.size());
}
++ii;
}
} else {
CHKERR skin.find_skin(0, els,
false, level_skin);
CHKERR pcomm->filter_pstatus(level_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
level_skin = subtract(level_skin, level0_skin);
adj, moab::Interface::UNION);
els = subtract(els, adj);
ele_to_refine.merge(els);
els,
SPACE_DIM - 1,
false, adj_edges, moab::Interface::UNION);
CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
}
CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr,
bit(
l),
true);
(boost::lexical_cast<std::string>(
l) +
"_ref_mesh.vtk").c_str(),
"VTK",
"");
(boost::lexical_cast<std::string>(
l) +
"_only_ref_mesh.vtk").c_str(),
"VTK", "");
};
auto mark_skins = [&](
auto l,
auto m) {
ents);
CHKERR skin.find_skin(0, ents,
false, level_skin);
CHKERR pcomm->filter_pstatus(level_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
level_skin = subtract(level_skin, level0_skin);
moab::Interface::UNION);
};
}
}
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p + 1; };
field_op_row->doWorkRhsHook = [](
DataOperator *op_ptr,
int side,
if (type == MBENTITYSET) {
<< "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
<<
" nb base functions integration points " << data.
getN().size1();
<< "\t" << CN::EntityTypeName(field_ent->getEntType());
}
}
};
}
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 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->divApproxVals = boost::make_shared<MatrixDouble>();
common_data_ptr->approxVals));
common_data_ptr->approxVals));
new OpErrorSkel<FIELD_DIM>(common_data_ptr));
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);
constexpr double eps = 1e-8;
"Not converged solution err = %6.4e", nrm2);
if (std::sqrt(array[0]) >
eps)
"Error in approximation err = %6.4e", std::sqrt(array[0]));
CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
}
auto rule = [](
int,
int,
int p) ->
int {
return -1; };
};
auto approx_vals = boost::make_shared<VectorDouble>();
double def_val[] = {0};
CHKERR moab.tag_get_handle(
"FIELD", 1, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
field_op_row->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
if (type == MBVERTEX) {
auto op_ptr =
base_op_ptr);
auto t_field = getFTensor0FromVec(*approx_vals);
auto nb_gauss_pts = op_ptr->getGaussPts().size2();
if (nb_gauss_pts != 3)
"Should be three guass pts.");
auto conn = op_ptr->getConn();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
const double v = t_field;
CHKERR moab.tag_set_data(th, &conn[gg], 1, &
v);
++t_field;
}
}
};
}
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
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.
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.
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.
constexpr double t
plate stiffness
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< MatrixDouble > divApproxVals
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 readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
static ApproxFieldFunctionDerivative< FIELD_DIM > divApproxFunction
MoFEM::Interface & mField
MoFEMErrorCode printResults()
[Check results]
MoFEMErrorCode runProblem()
[Run programme]
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.
FieldApproximationBase & getBase()
Get approximation base.
const VectorFieldEntities & getFieldEntities() const
get field entities
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FieldSpace & getSpace()
Get field space.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.