Implementation DG upwind method for advection/level set problem.
static char help[] =
"...\n\n";
IntegrationType::GAUSS;
constexpr bool debug =
false;
private:
struct SideData {
std::array<VectorInt, 2>
std::array<VectorInt, 2>
};
template <int SPACE_DIM>
static double get_level_set(
const double x,
const double y,
const double z);
boost::shared_ptr<FaceSideEle>
getSideFE(boost::shared_ptr<SideData> side_data_ptr);
boost::function<double(double, double, double)> level_fun =
boost::function<double(double, double, double)> vel_fun =
get_velocity_potential<SPACE_DIM>);
struct WrapperClass {
};
struct WrapperClassInitalSolution : public WrapperClass {
WrapperClassInitalSolution(boost::shared_ptr<double> max_ptr)
};
}
->synchroniseEntities(level);
}
return 0.05 * (*maxPtr);
}
private:
boost::shared_ptr<double>
maxPtr;
};
struct WrapperClassErrorProjection : public WrapperClass {
WrapperClassErrorProjection(boost::shared_ptr<double> max_ptr)
->synchroniseEntities(level);
}
double getThreshold(
const double max) {
return 0.05 * (*maxPtr); }
private:
boost::shared_ptr<double>
maxPtr;
};
struct OpRhsDomain;
struct OpLhsDomain;
struct OpRhsSkeleton;
I>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
I>::OpBaseTimesScalar<1>;
private:
boost::shared_ptr<double>
maxPtr;
};
template <>
double LevelSet::get_velocity_potential<2>(double x, double y, double z) {
return (x * x - 0.25) * (y * y - 0.25);
}
constexpr double xc = 0.1;
constexpr double yc = 0.;
constexpr double zc = 0.;
constexpr double r = 0.2;
return std::sqrt(pow(x - xc, 2) + pow(y - yc, 2) + pow(z - zc, 2)) -
r;
}
}
maxPtr = boost::make_shared<double>(0);
}
boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<VectorDouble> l_dot_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr);
private:
boost::shared_ptr<VectorDouble>
lPtr;
boost::shared_ptr<VectorDouble>
lDotPtr;
boost::shared_ptr<MatrixDouble>
velPtr;
};
boost::shared_ptr<MatrixDouble> vel_ptr);
private:
boost::shared_ptr<MatrixDouble>
velPtr;
};
boost::shared_ptr<FaceSideEle> side_fe_ptr);
private:
boost::shared_ptr<FaceSideEle>
};
boost::shared_ptr<FaceSideEle> side_fe_ptr);
private:
boost::shared_ptr<FaceSideEle>
};
int main(
int argc,
char *argv[]) {
try {
moab::Core moab_core;
moab::Interface &moab = moab_core;
DMType dm_name = "DMMOFEM";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "LevelSet"));
LogManager::setLog("LevelSet");
CHKERR level_set.runProblem();
}
return 0;
}
simple->getAddSkeletonFE() =
true;
simple->getAddBoundaryFE() =
true;
auto set_problem_bit = [&]() {
start_mask[s] = true;
#ifndef NDEBUG
CHKERR bit_mng->writeBitLevelByDim(
(proc_str + "level_base.vtk").c_str(), "VTK", "");
}
#endif
};
}
}
boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<VectorDouble> l_dot_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr)
lPtr(l_ptr), lDotPtr(l_dot_ptr), velPtr(vel_ptr) {}
boost::shared_ptr<MatrixDouble> vel_ptr)
velPtr(vel_ptr) {
this->sYmm = false;
}
boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr)
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr)
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
const auto nb_int_points = getGaussPts().size2();
const auto nb_base_func = data.
getN().size2();
auto t_l = getFTensor0FromVec(*
lPtr);
auto t_l_dot = getFTensor0FromVec(*
lDotPtr);
auto t_vel = getFTensor1FromMat<SPACE_DIM>(*
velPtr);
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != nb_int_points; ++gg) {
const auto alpha = t_w * getMeasure();
auto res0 =
alpha * t_l_dot;
t_res1(
i) = (
alpha * t_l) * t_vel(
i);
++t_w;
++t_l;
++t_l_dot;
++t_vel;
int rr = 0;
for (; rr != nb_dofs; ++rr) {
nf(rr) += res0 * t_base - t_res1(
i) * t_diff_base(
i);
++t_base;
++t_diff_base;
}
for (; rr < nb_base_func; ++rr) {
++t_base;
++t_diff_base;
}
}
}
const auto nb_int_points = getGaussPts().size2();
const auto nb_base_func = row_data.
getN().size2();
const auto nb_row_dofs = row_data.
getIndices().size();
const auto nb_col_dofs = col_data.
getIndices().size();
auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != nb_int_points; ++gg) {
const auto alpha = t_w * getMeasure();
const auto beta =
alpha * getTSa();
++t_w;
auto &mat = this->locMat;
int rr = 0;
for (; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; ++cc) {
mat(rr, cc) +=
(beta * t_row_base -
alpha * (t_row_diff_base(
i) * t_vel(
i))) *
t_col_base;
++t_col_base;
}
++t_row_base;
++t_row_diff_base;
}
for (; rr < nb_base_func; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_vel;
}
}
CHKERR loopSideFaces(
"dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop;
auto not_side = [](auto s) {
};
&*base_mat.data().begin());
};
if (in_the_loop > 0) {
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
if (nb_rows) {
resSkelton.resize(nb_rows, false);
resSkelton.clear();
const auto opposite_s0 = not_side(s0);
const auto sense_row = sideDataPtr->senseMap[s0];
#ifndef NDEBUG
const auto opposite_sense_row = sideDataPtr->senseMap[opposite_s0];
if (sense_row * opposite_sense_row > 0)
"Should be opposite sign");
#endif
const auto nb_row_base_functions =
sideDataPtr->rowBaseSideMap[s0].size2();
auto t_w = getFTensor0IntegrationWeight();
auto arr_t_l =
getFTensor0FromVec(sideDataPtr->lVec[
RIGHT_SIDE]));
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
auto next = [&]() {
for (auto &t_l : arr_t_l)
++t_l;
for (auto &t_vel : arr_t_vel)
++t_vel;
};
#ifndef NDEBUG
if (nb_gauss_pts != sideDataPtr->rowBaseSideMap[s0].size1())
"Inconsistent number of DOFs");
#endif
auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
const auto l_upwind = arr_t_l[l_upwind_side];
const auto res = t_w * dot * l_upwind;
next();
++t_w;
auto rr = 0;
for (; rr != nb_rows; ++rr) {
resSkelton[rr] += t_row_base * res;
++t_row_base;
}
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
}
}
CHKERR ::VecSetValues(getTSf(),
sideDataPtr->indicesRowSideMap[s0].size(),
&*sideDataPtr->indicesRowSideMap[s0].begin(),
&*resSkelton.begin(), ADD_VALUES);
}
}
}
}
CHKERR loopSideFaces(
"dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop;
auto not_side = [](auto s) {
};
&*base_mat.data().begin());
};
if (in_the_loop > 0) {
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
if (nb_rows) {
const auto opposite_s0 = not_side(s0);
const auto sense_row = sideDataPtr->senseMap[s0];
const auto nb_row_base_functions =
sideDataPtr->rowBaseSideMap[s0].size2();
const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
const auto sense_col = sideDataPtr->senseMap[s1];
matSkeleton.resize(nb_rows, nb_cols, false);
matSkeleton.clear();
auto t_w = getFTensor0IntegrationWeight();
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
auto next = [&]() {
for (auto &t_vel : arr_t_vel)
++t_vel;
};
auto t_row_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
const auto sense_upwind = sideDataPtr->senseMap[l_upwind_side];
auto res = t_w * dot;
next();
++t_w;
auto rr = 0;
if (s1 == l_upwind_side) {
for (; rr != nb_rows; ++rr) {
auto get_ntensor = [](
auto &base_mat,
auto gg,
auto bb) {
double *ptr = &base_mat(gg, bb);
};
auto t_col_base =
const auto res_row = res * t_row_base;
++t_row_base;
for (size_t cc = 0; cc != nb_cols; ++cc) {
matSkeleton(rr, cc) += res_row * t_col_base;
++t_col_base;
}
}
}
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
}
}
CHKERR ::MatSetValues(getTSB(),
sideDataPtr->indicesRowSideMap[s0].size(),
&*sideDataPtr->indicesRowSideMap[s0].begin(),
sideDataPtr->indicesColSideMap[s1].size(),
&*sideDataPtr->indicesColSideMap[s1].begin(),
&*matSkeleton.data().begin(), ADD_VALUES);
}
}
}
}
}
auto get_parent_vel_this = [&]() {
auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
parent_fe_ptr->getOpPtrVector().push_back(
"V", vel_ptr));
return parent_fe_ptr;
};
auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(
mField));
parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
}
return parents_elems_ptr_vec[0];
};
auto this_fe_ptr = get_parent_vel_this();
auto parent_fe_ptr = get_parents_vel_fe_ptr(this_fe_ptr);
}
pip->getOpDomainRhsPipeline().clear();
pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
pip->getDomainLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
pip->getDomainRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto l_ptr = boost::make_shared<VectorDouble>();
auto l_dot_ptr = boost::make_shared<VectorDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainRhsPipeline().push_back(
pip->getOpDomainRhsPipeline().push_back(
pip->getOpDomainRhsPipeline().push_back(
new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
}
boost::shared_ptr<FaceSideEle>
auto l_ptr = boost::make_shared<VectorDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
OpSideData(boost::shared_ptr<SideData> side_data_ptr)
sideDataPtr(side_data_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
sYmm = false;
}
if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
auto reset = [&](auto nb_in_loop) {
sideDataPtr->feSideHandle[nb_in_loop] = 0;
sideDataPtr->indicesRowSideMap[nb_in_loop].clear();
sideDataPtr->indicesColSideMap[nb_in_loop].clear();
sideDataPtr->rowBaseSideMap[nb_in_loop].clear();
sideDataPtr->colBaseSideMap[nb_in_loop].clear();
sideDataPtr->senseMap[nb_in_loop] = 0;
};
const auto nb_in_loop = getFEMethod()->nInTheLoop;
if (nb_in_loop == 0)
for (auto s : {0, 1})
reset(s);
sideDataPtr->currentFESide = nb_in_loop;
sideDataPtr->senseMap[nb_in_loop] = getSkeletonSense();
} else {
}
};
private:
boost::shared_ptr<SideData> sideDataPtr;
};
OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr)
sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
for (
auto t = moab::CN::TypeDimensionMap[
SPACE_DIM].first;
t <= moab::CN::TypeDimensionMap[
SPACE_DIM].second; ++
t)
sYmm = false;
}
if ((CN::Dimension(row_type) ==
SPACE_DIM) &&
const auto nb_in_loop = sideDataPtr->currentFESide;
sideDataPtr->feSideHandle[nb_in_loop] = getFEEntityHandle();
sideDataPtr->indicesRowSideMap[nb_in_loop] = row_data.
getIndices();
sideDataPtr->indicesColSideMap[nb_in_loop] = col_data.
getIndices();
sideDataPtr->rowBaseSideMap[nb_in_loop] = row_data.
getN();
sideDataPtr->colBaseSideMap[nb_in_loop] = col_data.
getN();
(sideDataPtr->lVec)[nb_in_loop] = *lPtr;
(sideDataPtr->velMat)[nb_in_loop] = *velPtr;
#ifndef NDEBUG
if ((sideDataPtr->lVec)[nb_in_loop].size() !=
(sideDataPtr->velMat)[nb_in_loop].size2())
"Wrong number of integaration pts %d != %d",
(sideDataPtr->lVec)[nb_in_loop].size(),
(sideDataPtr->velMat)[nb_in_loop].size2());
if ((sideDataPtr->velMat)[nb_in_loop].size1() !=
SPACE_DIM)
"Wrong size of velocity vector size = %d",
(sideDataPtr->velMat)[nb_in_loop].size1());
#endif
if (!nb_in_loop) {
(sideDataPtr->lVec)[1] = sideDataPtr->lVec[0];
(sideDataPtr->velMat)[1] = (sideDataPtr->velMat)[0];
} else {
#ifndef NDEBUG
if (sideDataPtr->rowBaseSideMap[0].size1() !=
sideDataPtr->rowBaseSideMap[1].size1()) {
"Wrong number of integration pt %d != %d",
sideDataPtr->rowBaseSideMap[0].size1(),
sideDataPtr->rowBaseSideMap[1].size1());
}
if (sideDataPtr->colBaseSideMap[0].size1() !=
sideDataPtr->colBaseSideMap[1].size1()) {
"Wrong number of integration pt");
}
#endif
}
} else {
}
};
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<VectorDouble> lPtr;
boost::shared_ptr<MatrixDouble> velPtr;
};
auto get_parent_this = [&]() {
auto parent_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
parent_fe_ptr->getOpPtrVector(), {potential_velocity_space, L2});
parent_fe_ptr->getOpPtrVector().push_back(
parent_fe_ptr->getOpPtrVector().push_back(
new OpSideDataOnParent(side_data_ptr, l_ptr, vel_ptr));
return parent_fe_ptr;
};
auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(
mField));
parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
}
return parents_elems_ptr_vec[0];
};
auto get_side_fe_ptr = [&]() {
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
auto this_fe_ptr = get_parent_this();
auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
side_fe_ptr->getOpPtrVector().push_back(new OpSideData(side_data_ptr));
side_fe_ptr->getOpPtrVector().push_back(
return side_fe_ptr;
};
return get_side_fe_ptr();
};
pip->getOpSkeletonRhsPipeline().clear();
pip->setSkeletonLhsIntegrationRule([](
int,
int,
int o) {
return 18; });
pip->setSkeletonRhsIntegrationRule([](
int,
int,
int o) {
return 18; });
pip->getSkeletonLhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
pip->getSkeletonRhsFE()->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto side_data_ptr = boost::make_shared<SideData>();
pip->getOpSkeletonRhsPipeline().push_back(
new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
pip->getOpSkeletonLhsPipeline().push_back(
}
OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
boost::shared_ptr<SideData> side_data_ptr,
sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
errorSumPtr(error_sum_ptr), thError(th_error) {}
CHKERR loopSideFaces(
"dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop;
auto not_side = [](auto s) {
};
auto nb_gauss_pts = getGaussPts().size2();
auto arr_t_l =
getFTensor0FromVec(sideDataPtr->lVec[
RIGHT_SIDE]));
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
auto next = [&]() {
for (auto &t_l : arr_t_l)
++t_l;
for (auto &t_vel : arr_t_vel)
++t_vel;
};
double e = 0;
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
e += t_w * getMeasure() *
next();
++t_w;
}
e = std::sqrt(e);
moab::Interface &moab =
getNumeredEntFiniteElementPtr()->getBasicDataPtr()->moab;
const void *tags_ptr[2];
CHKERR moab.tag_get_by_ptr(thError, sideDataPtr->feSideHandle.data(), 2,
tags_ptr);
for (auto ff : {0, 1}) {
*((double *)tags_ptr[ff]) += e;
}
CHKERR VecSetValue(errorSumPtr, 0, e, ADD_VALUES);
};
}
private:
boost::shared_ptr<FaceSideEle> sideFEPtr;
boost::shared_ptr<SideData> sideDataPtr;
Tag thError;
};
Tag th_error;
double def_val = 0;
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_val);
auto clear_tags = [&]() {
double zero;
};
auto evaluate_error = [&]() {
auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
skel_fe->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
auto side_data_ptr = boost::make_shared<SideData>();
skel_fe->getOpPtrVector().push_back(
new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
};
simple->getSkeletonFEName(), skel_fe);
};
auto assemble_and_sum = [](auto vec) {
double sum;
return sum;
};
auto propagate_error_to_parents = [&]() {
auto fe_ptr = boost::make_shared<FEMethod>();
fe_ptr->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
fe_ptr->preProcessHook = []() { return 0; };
fe_ptr->postProcessHook = []() { return 0; };
fe_ptr->operatorHook = [&]() {
auto fe_ent = fe_ptr->numeredEntFiniteElementPtr->getEnt();
auto parent = fe_ptr->numeredEntFiniteElementPtr->getParentEnt();
auto th_parent = fe_ptr->numeredEntFiniteElementPtr->getBasicDataPtr()
->th_RefParentHandle;
double error;
CHKERR moab.tag_get_data(th_error, &fe_ent, 1, &error);
[&](auto fe_ent, auto error) {
double *e_ptr;
CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
(const void **)&e_ptr);
(*e_ptr) += error;
CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
if (parent != fe_ent && parent)
CHKERR add_error(parent, *e_ptr);
};
CHKERR add_error(parent, error);
};
fe_ptr);
};
return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
}
DivergenceVol(boost::shared_ptr<VectorDouble> l_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr,
divVec(div_vec) {}
if (nb_dofs) {
const auto nb_gauss_pts = getGaussPts().size2();
const auto t_w = getFTensor0IntegrationWeight();
auto t_l = getFTensor0FromVec(*lPtr);
auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
double div = 0;
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
for (int rr = 0; rr != nb_dofs; ++rr) {
div += getMeasure() * t_w * t_l * (t_diff(
i) * t_vel(
i));
++t_diff;
}
++t_w;
++t_l;
++t_vel;
}
CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
}
}
private:
boost::shared_ptr<VectorDouble> lPtr;
boost::shared_ptr<MatrixDouble> velPtr;
};
DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr,
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
&*base_mat.data().begin());
};
auto not_side = [](auto s) {
};
CHKERR loopSideFaces(
"dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop;
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
const auto nb_dofs = sideDataPtr->indicesRowSideMap[s0].size();
if (nb_dofs) {
auto t_base =
get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
auto nb_row_base_functions = sideDataPtr->rowBaseSideMap[s0].size2();
auto side_sense = sideDataPtr->senseMap[s0];
auto opposite_s0 = not_side(s0);
auto arr_t_l =
getFTensor0FromVec(sideDataPtr->lVec[
RIGHT_SIDE]));
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
LEFT_SIDE]),
getFTensor1FromMat<SPACE_DIM>(sideDataPtr->velMat[
RIGHT_SIDE]));
auto next = [&]() {
for (auto &t_l : arr_t_l)
++t_l;
for (auto &t_vel : arr_t_vel)
++t_vel;
};
double div = 0;
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
const auto dot = (t_normal(
i) * t_vel(
i)) * side_sense;
const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
const auto l_upwind =
arr_t_l[l_upwind_side];
auto res = t_w * l_upwind * dot;
++t_w;
next();
int rr = 0;
for (; rr != nb_dofs; ++rr) {
div += t_base * res;
++t_base;
}
for (; rr < nb_row_base_functions; ++rr) {
++t_base;
}
}
CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
}
if (!in_the_loop)
break;
}
}
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<FaceSideEle> sideFEPtr;
boost::shared_ptr<MatrixDouble> velPtr;
};
auto vol_fe = boost::make_shared<DomainEle>(
mField);
auto skel_fe = boost::make_shared<BoundaryEle>(
mField);
vol_fe->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
skel_fe->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
auto l_ptr = boost::make_shared<VectorDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
auto side_data_ptr = boost::make_shared<SideData>();
vol_fe->getOpPtrVector(), {potential_velocity_space, L2});
vol_fe->getOpPtrVector().push_back(
vol_fe->getOpPtrVector().push_back(
new DivergenceVol(l_ptr, vel_ptr, div_vol_vec));
skel_fe->getOpPtrVector().push_back(
new DivergenceSkeleton(side_data_ptr, side_fe_ptr, div_skel_vec));
[](double x, double y, double) { return x - y; });
[](double x, double y, double) { return x - y; });
vol_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
skel_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto assemble_and_sum = [](auto vec) {
double sum;
return sum;
};
auto div_vol = assemble_and_sum(div_vol_vec);
auto div_skel = assemble_and_sum(div_skel_vec);
auto eps = std::abs((div_vol - div_skel) / (div_vol + div_skel));
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence volume: " << div_vol;
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Testing divergence skeleton: " << div_skel
<<
" relative difference: " <<
eps;
constexpr double eps_err = 1e-6;
"No consistency between skeleton integral and volume integral");
};
auto post_proc = [&](auto dm, auto f_res, auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
auto l_vec = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{},
{}, {})
);
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
constexpr double eps = 1e-4;
auto x =
opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}, {
"V", {-1, 1}}});
auto dot_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
auto diff_x = opt->setRandomFields(
simple->getDM(), {{
"L", {-1, 1}}});
auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
auto rhs_pipeline) {
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
CHKERR post_proc(
simple->getDM(), diff_res,
"tangent_op_error.h5m");
}
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
"Test consistency of tangent matrix %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err)
"Norm of directional derivative too large err = %3.4e", fnorm);
};
CHKERR test_domain_ops(
simple->getDomainFEName(), pip->getDomainLhsFE(),
pip->getDomainRhsFE());
CHKERR test_domain_ops(
simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
pip->getSkeletonRhsFE());
};
boost::function<double(double, double, double)> level_fun) {
boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
auto swap_fe = [&]() {
lhs_fe.swap(pip->getDomainLhsFE());
rhs_fe.swap(pip->getDomainRhsFE());
};
swap_fe();
pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
remove_mask.flip();
CHKERR prb_mng->removeDofsOnEntities(
"LEVELSET_POJECTION",
"L",
auto test_bit_ele = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainLhsPipeline().push_back(
new OpMassLL(
"L",
"L"));
pip->getOpDomainRhsPipeline().push_back(
new OpSourceL(
"L", level_fun));
auto ksp = pip->createKSP(sub_dm);
CHKERR KSPSetFromOptions(ksp);
CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
#ifndef NDEBUG
auto fe_meshset =
std::vector<Tag> tags{th_error};
"PARALLEL=WRITE_PART", &fe_meshset, 1,
&*tags.begin(), tags.size());
#endif
auto post_proc = [&](auto dm, auto out_name, auto th_error) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
post_proc_fe->setTagsToTransfer({th_error});
post_proc_fe->exeTestHook = test_bit_ele;
auto l_vec = boost::make_shared<VectorDouble>();
auto l_grad_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{{"GradL", l_grad_mat}},
{}, {})
);
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
CHKERR post_proc(sub_dm,
"initial_level_set.h5m", th_error);
swap_fe();
}
boost::function<double(double, double, double)> vel_fun) {
boost::shared_ptr<FEMethod> lhs_fe = boost::make_shared<DomainEle>(
mField);
boost::shared_ptr<FEMethod> rhs_fe = boost::make_shared<DomainEle>(
mField);
auto swap_fe = [&]() {
lhs_fe.swap(pip->getDomainLhsFE());
rhs_fe.swap(pip->getDomainRhsFE());
};
swap_fe();
pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
pip->setDomainRhsIntegrationRule([](
int,
int,
int o) {
return 3 *
o; });
remove_mask.flip();
CHKERR prb_mng->removeDofsOnEntities(
"VELOCITY_PROJECTION",
"V",
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
};
pip->getDomainLhsFE()->exeTestHook = test_bit;
pip->getDomainRhsFE()->exeTestHook = test_bit;
pip->getOpDomainLhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainRhsPipeline(), {potential_velocity_space, L2});
pip->getOpDomainLhsPipeline().push_back(
new OpMassVV(
"V",
"V"));
pip->getOpDomainRhsPipeline().push_back(
new OpSourceV(
"V", vel_fun));
auto ksp = pip->createKSP(sub_dm);
CHKERR KSPSetFromOptions(ksp);
CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
auto post_proc = [&](auto dm, auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
post_proc_fe->exeTestHook = test_bit;
post_proc_fe->getOpPtrVector(), {potential_velocity_space});
auto l_vec = boost::make_shared<VectorDouble>();
auto potential_vec = boost::make_shared<VectorDouble>();
auto velocity_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"VelocityPotential", potential_vec}},
{{"Velocity", velocity_mat}},
{}, {})
);
} else {
"3d case not implemented");
}
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
CHKERR post_proc(sub_dm,
"initial_velocity_potential.h5m");
swap_fe();
}
remove_mask.flip();
remove_mask);
auto add_post_proc_fe = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
Tag th_error;
double def_val = 0;
"Error", 1, MB_TYPE_DOUBLE, th_error, MB_TAG_CREAT | MB_TAG_SPARSE,
&def_val);
post_proc_fe->setTagsToTransfer({th_error});
post_proc_fe->exeTestHook = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto l_vec = boost::make_shared<VectorDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector(), {H1});
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{{"V", vel_ptr}},
{}, {})
);
return post_proc_fe;
};
auto post_proc_fe = add_post_proc_fe();
auto set_time_monitor = [&](auto dm, auto ts) {
auto monitor_ptr = boost::make_shared<FEMethod>();
monitor_ptr->preProcessHook = []() { return 0; };
monitor_ptr->operatorHook = []() { return 0; };
monitor_ptr->postProcessHook = [&]() {
if (!post_proc_fe)
"Null pointer for post proc element");
post_proc_fe);
CHKERR post_proc_fe->writeFile(
"level_set_" +
boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
};
boost::shared_ptr<FEMethod> null;
null, null);
return monitor_ptr;
};
auto ts = pip->createTSIM(sub_dm);
auto set_solution = [&](auto ts) {
};
auto monitor_pt = set_time_monitor(sub_dm, ts);
CHKERR TSSetIJacobian(ts,
B,
B, TsSetIJacobian,
nullptr);
auto ts_pre_step = [](TS ts) {
MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
auto set_solution = [&](auto ts) {
DM dm;
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
<< "Problem " << prb_ptr->getName() << " solution vector norm "
<< get_norm(x);
};
auto refine_and_project = [&](auto ts) {
DM dm;
remove_mask
.flip();
->removeDofsOnEntities(
"ADVECTION",
"L",
BitRefLevel().set(),
remove_mask);
};
auto ts_reset_theta = [&](auto ts) {
DM dm;
CHKERR TSSetIJacobian(ts,
B,
B, TsSetIJacobian,
nullptr);
};
CHKERR refine_and_project(ts);
};
auto ts_post_step = [](TS ts) { return 0; };
CHKERR TSSetPreStep(ts, ts_pre_step);
CHKERR TSSetPostStep(ts, ts_post_step);
}
ParallelComm *pcomm =
meshset_ptr->get_ptr(), 1);
};
auto get_refined_elements_meshset = [&](
auto bit,
auto mask) {
Tag th_error;
"get error handle");
std::vector<double> errors(fe_ents.size());
"get tag data");
auto it = std::max_element(errors.begin(), errors.end());
double max;
MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Max error: " << max;
auto threshold = wp.getThreshold(max);
std::vector<EntityHandle> fe_to_refine;
fe_to_refine.reserve(fe_ents.size());
auto fe_it = fe_ents.begin();
auto error_it = errors.begin();
for (
auto i = 0;
i != fe_ents.size(); ++
i) {
if (*error_it > threshold) {
fe_to_refine.push_back(*fe_it);
}
++fe_it;
++error_it;
}
ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
auto get_neighbours_by_bridge_vertices = [&](auto &&ents) {
moab::Interface::UNION);
return ents;
};
ents = get_neighbours_by_bridge_vertices(ents);
#ifndef NDEBUG
"add entities to meshset");
(proc_str + "_fe_to_refine.vtk").c_str(), "VTK", "",
meshset_ptr->get_ptr(), 1);
}
#endif
return ents;
};
auto refine_mesh = [&](
auto l,
auto &&fe_to_refine) {
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
fe_to_refine = intersect(level_ents, fe_to_refine);
level_ents = subtract(level_ents, fe_to_refine);
Range fe_to_refine_children;
bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
fe_to_refine_children =
fe_to_refine_children.subset_by_dimension(
SPACE_DIM);
level_ents.merge(fe_to_refine_children);
auto fix_neighbour_level = [&](auto ll) {
auto level_ll = level_ents;
level_ll);
CHKERR skin.find_skin(0, level_ll,
false, skin_edges);
for (auto lll = 0; lll <= ll; ++lll) {
CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
}
for (auto lll = 0; lll <= ll - 2; ++lll) {
bad_bit[lll] = true;
}
skin_adj_ents,
moab::Interface::UNION);
skin_adj_ents);
skin_adj_ents = intersect(skin_adj_ents, level_ents);
if (!skin_adj_ents.empty()) {
level_ents = subtract(level_ents, skin_adj_ents);
Range skin_adj_ents_children;
bit_mng->updateRangeByChildren(skin_adj_ents, skin_adj_ents_children);
level_ents.merge(skin_adj_ents_children);
}
};
level_ents);
if (d == 0) {
level_ents.subset_by_dimension(
SPACE_DIM), level_ents,
true);
} else {
level_ents.subset_by_dimension(
SPACE_DIM), d,
false, level_ents,
moab::Interface::UNION);
}
}
level_ents);
#ifndef NDEBUG
CHKERR bit_mng->writeBitLevelByDim(
(boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_ref_mesh.vtk")
.c_str(),
"VTK", "");
#endif
};
auto set_skelton_bit = [&](
auto l) {
Range level_edges_parents;
CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
level_edges_parents =
level_edges_parents.subset_by_dimension(
SPACE_DIM - 1);
CHKERR bit_mng->filterEntitiesByRefLevel(
auto parent_skeleton = intersect(level_edges, level_edges_parents);
auto skeleton = subtract(level_edges, level_edges_parents);
moab::Interface::UNION);
#ifndef NDEBUG
CHKERR bit_mng->writeBitLevel(
(boost::lexical_cast<std::string>(
l) +
"_" + proc_str +
"_skeleton.vtk")
.c_str(),
"VTK", "");
#endif
};
start_mask[s] = true;
CHKERR bit_mng->lambdaBitRefLevel(
true);
level0);
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Process level: " <<
l;
CHKERR refine_mesh(
l + 1, get_refined_elements_meshset(
CHKERR wp.setAggregateBit(*
this,
l + 1);
}
}
auto lhs_fe = boost::make_shared<DomainEle>(
mField);
auto rhs_fe_prj = boost::make_shared<DomainEle>(
mField);
auto rhs_fe_current = boost::make_shared<DomainEle>(
mField);
lhs_fe->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
rhs_fe_prj->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
rhs_fe_current->getRuleHook = [](
int,
int,
int o) {
return 3 *
o; };
current_ents);
prj_ents);
CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
}
current_ents = subtract(current_ents, prj_ents);
auto test_mesh_bit = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto test_prj_bit = [&](
FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
};
auto test_current_bit = [&](
FEMethod *fe_ptr) {
return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
};
lhs_fe->exeTestHook = test_mesh_bit;
rhs_fe_prj->exeTestHook = test_prj_bit;
rhs_fe_current->exeTestHook = test_current_bit;
remove_mask.flip();
CHKERR prb_mng->removeDofsOnEntities(
"DG_PROJECTION",
"L",
BitRefLevel().set(), remove_mask,
nullptr, 0,
lhs_fe->getOpPtrVector(), {potential_velocity_space, L2});
lhs_fe->getOpPtrVector().push_back(
new OpMassLL(
"L",
"L"));
auto l_vec = boost::make_shared<VectorDouble>();
auto set_prj_from_child = [&](auto rhs_fe_prj) {
rhs_fe_prj->getOpPtrVector(), {potential_velocity_space, L2});
rhs_fe_prj->getOpPtrVector().push_back(
auto get_parent_this = [&]() {
auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
fe_parent_this->getOpPtrVector().push_back(
return fe_parent_this;
};
auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(
mField));
parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
}
return parents_elems_ptr_vec[0];
};
auto this_fe_ptr = get_parent_this();
auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
rhs_fe_prj->getOpPtrVector().push_back(
};
auto set_prj_from_parent = [&](auto rhs_fe_current) {
auto get_parent_this = [&]() {
auto fe_parent_this = boost::make_shared<DomianParentEle>(
mField);
fe_parent_this->getOpPtrVector().push_back(
return fe_parent_this;
};
auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(
mField));
parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
}
return parents_elems_ptr_vec[0];
};
auto this_fe_ptr = get_parent_this();
auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
false);
l_vec->clear();
return 0;
};
rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
rhs_fe_current->getOpPtrVector().push_back(
rhs_fe_current->getOpPtrVector().push_back(
new OpScalarFieldL(
"L", l_vec));
};
set_prj_from_child(rhs_fe_prj);
set_prj_from_parent(rhs_fe_current);
boost::shared_ptr<FEMethod> null_fe;
lhs_fe, null_fe, null_fe);
null_fe, null_fe);
rhs_fe_current, null_fe, null_fe);
CHKERR KSPSetFromOptions(ksp);
CHKERR VecGhostUpdateBegin(
L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
L, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG(
"LevelSet", Sev::inform) <<
"Error indicator " << error;
auto post_proc = [&](auto dm, auto out_name, auto th_error) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
post_proc_fe->setTagsToTransfer({th_error});
post_proc_fe->exeTestHook = test_mesh_bit;
auto l_vec = boost::make_shared<VectorDouble>();
auto l_grad_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector(), {potential_velocity_space, L2});
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{{"GradL", l_grad_mat}},
{}, {})
);
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
CHKERR post_proc(sub_dm,
"dg_projection.h5m", th_error);
}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ 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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto get_ntensor(T &base_mat)
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
constexpr int skeleton_bit
skeleton elemets bit
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
constexpr int current_bit
dofs bit used to do calculations
constexpr FieldSpace potential_velocity_space
FTensor::Index< 'i', SPACE_DIM > i
constexpr int SPACE_DIM
[Define dimension]
FaceSideEle::UserDataOperator FaceSideEleOp
constexpr int aggregate_projection_bit
all bits for projection problem
constexpr IntegrationType I
DomainEle::UserDataOperator DomainEleOp
LevelSet * level_set_raw_ptr
constexpr int aggregate_bit
all bits for advection problem
constexpr int projection_bit
constexpr size_t potential_velocity_field_dim
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
constexpr auto make_array(Arg &&...arg)
Create Array.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
constexpr double t
plate stiffness
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
OpLhsDomain(const std::string field_name, boost::shared_ptr< MatrixDouble > vel_ptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > velPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
boost::shared_ptr< SideData > sideDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpLhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
boost::shared_ptr< VectorDouble > lPtr
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > velPtr
OpRhsDomain(const std::string field_name, boost::shared_ptr< VectorDouble > l_ptr, boost::shared_ptr< VectorDouble > l_dot_ptr, boost::shared_ptr< MatrixDouble > vel_ptr)
boost::shared_ptr< VectorDouble > lDotPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< SideData > sideDataPtr
OpRhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
std::array< EntityHandle, 2 > feSideHandle
int currentFESide
current side counter
std::array< MatrixDouble, 2 > colBaseSideMap
std::array< MatrixDouble, 2 > rowBaseSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
std::array< int, 2 > senseMap
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
MoFEMErrorCode setBits(LevelSet &level_set, int l)
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
boost::shared_ptr< double > maxPtr
virtual double getThreshold(const double max)=0
virtual MoFEMErrorCode runCalcs(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)=0
virtual MoFEMErrorCode setBits(LevelSet &level_set, int l)=0
double getThreshold(const double max)
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l)
boost::shared_ptr< double > maxPtr
MoFEMErrorCode setBits(LevelSet &level_set, int l)
MoFEMErrorCode runCalcs(LevelSet &level_set, int l)
MoFEM::Interface & mField
integrate skeleton operators on khs
MoFEMErrorCode solveAdvection()
solve advection problem
boost::shared_ptr< double > maxPtr
std::array< MatrixDouble, 2 > MatSideArray
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< 1, 1 > OpMassLL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
MoFEMErrorCode readMesh()
read mesh
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
MoFEMErrorCode setupProblem()
create fields, and set approximation order
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
MoFEMErrorCode runProblem()
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< 1 > OpScalarFieldL
static double get_velocity_potential(double x, double y, double z)
advection velocity field
std::array< VectorDouble, 2 > VecSideArray
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< SPACE_DIM >)
initialise potential velocity field
static double get_level_set(const double x, const double y, const double z)
inital level set, i.e. advected filed
MoFEMErrorCode refineMesh(WrapperClass &&wp)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpSource< 1, 1 > OpSourceL
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
std::tuple< double, Tag > evaluateError()
evaluate error
MoFEMErrorCode testSideFE()
test integration side elements
Add operators pushing bases from local to physical configuration.
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::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
Calculate curl of vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get rate of scalar field at integration points.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
intrusive_ptr for managing petsc objects
Interface for Time Stepping (TS) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator the left hand side matrix.