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,
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));
if (op == DomainEleOp::OPSPACE) {
fe_ptr_current->getOpPtrVector(), {H1});
}
add_parent_level(
boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
fe_ptr_current),
level - 1);
if (op == DomainEleOp::OPSPACE) {
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:
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> divApproxVals;
};
template <
int FIELD_DIM>
struct OpError;
template <int FIELD_DIM> struct OpErrorSkel;
};
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
if (
const size_t nb_dofs = data.
getIndices().size()) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad_val =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
auto t_coords = getFTensor1CoordsAtGaussPts();
nf.clear();
const double volume = getMeasure();
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;
}
const int index = 0;
CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
}
}
};
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_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 =
Skinner skin(&mField.get_moab());
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Dim " << simpleInterface->getDim();
auto &moab = mField.get_moab();
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);
CHKERR mField.get_moab().get_adjacencies(
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);
CHKERR mField.get_moab().get_adjacencies(level_skin, 0,
false, level_skin,
moab::Interface::UNION);
};
}
simpleInterface->getBitRefLevel() = bit_sum;
simpleInterface->getBitRefLevelMask() =
BitRefLevel().set();
}
simpleInterface->getParentAdjacencies() = true;
CHKERR simpleInterface->setUp();
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 * p + 1; };
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainLhsFE(),
DomainEleOp::OPSPACE,
QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainLhsFE(),
DomainEleOp::OPROW,
QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainLhsFE(),
DomainEleOp::OPCOL,
QUIET, Sev::noisy);
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());
}
}
};
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPSPACE,
NOISY, Sev::verbose);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPROW,
NOISY, Sev::noisy);
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
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>();
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPSPACE,
QUIET, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPROW,
VERBOSE, Sev::noisy);
common_data_ptr->approxVals));
BoundaryEleOp::OPSPACE,
QUIET, Sev::noisy);
BoundaryEleOp::OPROW,
VERBOSE, Sev::noisy);
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>();
auto &moab = mField.get_moab();
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,
auto op_ptr =
base_op_ptr);
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;
}
}
};
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPSPACE,
VERBOSE, Sev::noisy);
set_parent_dofs<DomainParentEle>(mField, pipeline_mng->
getDomainRhsFE(),
DomainEleOp::OPROW,
VERBOSE, Sev::noisy);
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
}
}