v0.14.0
level_set.cpp

Implementation DG upwind method for advection/level set problem

Date
2022-12-15
Copyright
Copyright (c) 2022
/**
* @file level_set.cpp
* @example level_set.cpp
* @brief Implementation DG upwind method for advection/level set problem
* @date 2022-12-15
*
* @copyright Copyright (c) 2022
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
//! [Define dimension]
// We can have 2D elements embedded in 3D space or 2D element embedded in 2D
constexpr int FE_DIM = EXECUTABLE_DIMENSION; //< dimension of the element
constexpr int SPACE_DIM = FE_DIM; //< dimension of the space
constexpr int DIM1 = 1;
constexpr int DIM2 = 1;
// \todo We do not have implemented embedding of the edge elements in 3D space,
// thus integration on skeleton will not work for
// case FE_DIM = 2 and SPACE_DIM = 3. Since FE_DIM = 2, skeleton is 1D, i.e.
// edge elements.
FTensor::Index<'I', DIM1> I;
FTensor::Index<'J', DIM1> J;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType G =
IntegrationType::GAUSS; //< selected integration type
using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<FE_DIM>::DomainEle;
using DomianParentEle =
PipelineManager::ElementsAndOpsByDim<FE_DIM>::DomianParentEle;
using BoundaryEle = PipelineManager::ElementsAndOpsByDim<FE_DIM>::BoundaryEle;
using FaceSideEle = PipelineManager::ElementsAndOpsByDim<FE_DIM>::FaceSideEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using FaceSideEleOp = FaceSideEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
constexpr FieldSpace potential_velocity_space = FE_DIM == 2 ? H1 : HCURL;
constexpr size_t potential_velocity_field_dim = FE_DIM == 2 ? 1 : 3;
// #ifndef NDEBUG
constexpr bool debug = true;
// #else
// constexpr bool debug = true;
// #endif
constexpr int nb_levels = 3; //< number of refinement levels
constexpr int start_bit =
nb_levels + 1; //< first refinement level for computational mesh
constexpr int current_bit =
2 * start_bit + 1; ///< dofs bit used to do calculations
constexpr int skeleton_bit = 2 * start_bit + 2; ///< skeleton elements bit
constexpr int aggregate_bit =
2 * start_bit + 3; ///< all bits for advection problem
constexpr int projection_bit =
2 * start_bit + 4; //< bit from which data are projected
constexpr int aggregate_projection_bit =
2 * start_bit + 5; ///< all bits for projection problem
struct LevelSet {
LevelSet(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
using MatSideArray = std::array<MatrixDouble, 2>;
/**
* @brief data structure carrying information on skeleton on both sides.
*
*/
struct SideData {
// data for skeleton computation
std::array<EntityHandle, 2> feSideHandle;
std::array<VectorInt, 2>
indicesRowSideMap; ///< indices on rows for left hand-side
std::array<VectorInt, 2>
indicesColSideMap; ///< indices on columns for left hand-side
std::array<MatrixDouble, 2> rowBaseSideMap; // base functions on rows
std::array<MatrixDouble, 2> colBaseSideMap; // base function on columns
std::array<int, 2> senseMap; // orientation of local element edge/face in
// respect to global orientation of edge/face
MatSideArray lVec; //< Values of level set field
MatSideArray velMat; //< Values of velocity field
int currentFESide; ///< current side counter
};
/**
* @brief advection velocity field
*
* \note in current implementation is assumed that advection field has zero
* normal component.
*
* \note function define a vector velocity potential field, curl of potential
* field gives velocity, thus velocity is divergence free.
*
* @tparam FE_DIM
* @param x
* @param y
* @param z
* @return auto
*/
template <int FE_DIM>
static double get_velocity_potential(double x, double y, double z);
/**
* @brief inital level set, i.e. advected filed
*
* @param x
* @param y
* @param z
* @return double
*/
static double get_level_set(const double x, const double y, const double z);
/**
* @brief read mesh
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode readMesh();
/**
* @brief create fields, and set approximation order
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode setUpProblem();
/**
* @brief push operators to integrate operators on domain
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode pushOpDomain();
/**
* @brief evaluate error
*
* @return MoFEMErrorCode
*/
std::tuple<double, Tag> evaluateError();
/**
* @brief Get operator calculating velocity on coarse mesh
*
* @param vel_ptr
* @return DomainEleOp*
*/
ForcesAndSourcesCore::UserDataOperator *
getZeroLevelVelOp(boost::shared_ptr<MatrixDouble> vel_ptr);
/**
* @brief create side element to assemble data from sides
*
* @param side_data_ptr
* @return boost::shared_ptr<FaceSideEle>
*/
boost::shared_ptr<FaceSideEle>
getSideFE(boost::shared_ptr<SideData> side_data_ptr);
/**
* @brief push operator to integrate on skeleton
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode pushOpSkeleton();
/**
* @brief test integration side elements
*
* Check consistency between volume and skeleton integral.
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode testSideFE();
/**
* @brief test consistency between tangent matrix and the right hand side
* vectors
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode testOp();
/**
* @brief initialise field set
*
* @param level_fun
* @return MoFEMErrorCode
*/
MoFEMErrorCode initialiseFieldLevelSet(
boost::function<double(double, double, double)> level_fun =
get_level_set);
/**
* @brief initialise potential velocity field
*
* @param vel_fun
* @return MoFEMErrorCode
*/
MoFEMErrorCode initialiseFieldVelocity(
boost::function<double(double, double, double)> vel_fun =
get_velocity_potential<FE_DIM>);
/**
* @brief dg level set projection
*
* @param prj_bit
* @param mesh_bit
* @return MoFEMErrorCode
*/
MoFEMErrorCode dgProjection(const int prj_bit = projection_bit);
/**
* @brief solve advection problem
*
* @return * MoFEMErrorCode
*/
MoFEMErrorCode solveAdvection();
/**
* @brief Wrapper executing stages while mesh refinement
*/
struct WrapperClass {
WrapperClass() = default;
/**
* @brief Set bit ref level to problem
*/
virtual MoFEMErrorCode setBits(LevelSet &level_set, int l) = 0;
/**
* @brief Run calculations
*/
virtual MoFEMErrorCode runCalcs(LevelSet &level_set, int l) = 0;
/**
* @brief Add bit to current element, so it aggregate all previious current
* elements
*/
virtual MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l) = 0;
virtual double getThreshold(const double max) = 0;
};
/**
* @brief Used to execute inital mesh approximation while mesh refinement
*
*/
struct WrapperClassInitalSolution : public WrapperClass {
WrapperClassInitalSolution(boost::shared_ptr<double> max_ptr)
: WrapperClass(), maxPtr(max_ptr) {}
MoFEMErrorCode setBits(LevelSet &level_set, int l) {
MoFEMFunctionBegin;
auto simple = level_set.mField.getInterface<Simple>();
simple->getBitRefLevel() =
BitRefLevel().set(skeleton_bit) | BitRefLevel().set(aggregate_bit);
simple->getBitRefLevelMask() = BitRefLevel().set();
simple->reSetUp(true);
MoFEMFunctionReturn(0);
};
MoFEMErrorCode runCalcs(LevelSet &level_set, int l) {
MoFEMFunctionBegin;
CHKERR level_set.initialiseFieldLevelSet();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l) {
auto bit_mng = level_set.mField.getInterface<BitRefManager>();
auto set_bit = [](auto l) { return BitRefLevel().set(l); };
MoFEMFunctionBegin;
Range level;
CHKERR bit_mng->getEntitiesByRefLevel(set_bit(start_bit + l),
BitRefLevel().set(), level);
CHKERR level_set.mField.getInterface<CommInterface>()
->synchroniseEntities(level);
CHKERR bit_mng->setNthBitRefLevel(current_bit, false);
CHKERR bit_mng->setNthBitRefLevel(level, current_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level, aggregate_bit, true);
MoFEMFunctionReturn(0);
}
double getThreshold(const double max) {
*maxPtr = std::max(*maxPtr, max);
return 0.05 * (*maxPtr);
}
private:
boost::shared_ptr<double> maxPtr;
};
/**
* @brief Use peculated errors on all levels while mesh projection
*
*/
struct WrapperClassErrorProjection : public WrapperClass {
WrapperClassErrorProjection(boost::shared_ptr<double> max_ptr)
: maxPtr(max_ptr) {}
MoFEMErrorCode setBits(LevelSet &level_set, int l) { return 0; };
MoFEMErrorCode runCalcs(LevelSet &level_set, int l) { return 0; }
MoFEMErrorCode setAggregateBit(LevelSet &level_set, int l) {
auto bit_mng = level_set.mField.getInterface<BitRefManager>();
auto set_bit = [](auto l) { return BitRefLevel().set(l); };
MoFEMFunctionBegin;
Range level;
CHKERR bit_mng->getEntitiesByRefLevel(set_bit(start_bit + l),
BitRefLevel().set(), level);
CHKERR level_set.mField.getInterface<CommInterface>()
->synchroniseEntities(level);
CHKERR bit_mng->setNthBitRefLevel(current_bit, false);
CHKERR bit_mng->setNthBitRefLevel(level, current_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level, aggregate_bit, true);
MoFEMFunctionReturn(0);
}
double getThreshold(const double max) { return 0.05 * (*maxPtr); }
private:
boost::shared_ptr<double> maxPtr;
};
MoFEMErrorCode refineMesh(WrapperClass &&wp);
struct OpRhsDomain; ///< integrate volume operators on rhs
struct OpLhsDomain; ///< integrate volume operator on lhs
struct OpRhsSkeleton; ///< integrate skeleton operators on rhs
struct OpLhsSkeleton; ///< integrate skeleton operators on khs
// Main interfaces
MoFEM::Interface &mField;
using AssemblyDomainEleOp =
FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
using OpMassLL = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
G>::OpMass<1, DIM1 * DIM2>;
using OpSourceL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
G>::OpSource<1, DIM1 * DIM2>;
using OpMassVV = FormsIntegrators<DomainEleOp>::Assembly<A>::BiLinearForm<
G>::OpMass<potential_velocity_field_dim, potential_velocity_field_dim>;
using OpSourceV = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
G>::OpSource<potential_velocity_field_dim, potential_velocity_field_dim>;
using OpScalarFieldL = FormsIntegrators<DomainEleOp>::Assembly<A>::LinearForm<
G>::OpBaseTimesVector<1, DIM1 * DIM2, 1>;
using AssemblyBoundaryEleOp =
FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
enum ElementSide { LEFT_SIDE = 0, RIGHT_SIDE = 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);
}
double LevelSet::get_level_set(const double x, const double y, const double z) {
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;
}
MoFEMErrorCode LevelSet::runProblem() {
MoFEMFunctionBegin;
CHKERR readMesh();
CHKERR setUpProblem();
if constexpr (debug) {
CHKERR testSideFE();
CHKERR testOp();
}
CHKERR initialiseFieldVelocity();
maxPtr = boost::make_shared<double>(0);
CHKERR refineMesh(WrapperClassInitalSolution(maxPtr));
auto simple = mField.getInterface<Simple>();
simple->getBitRefLevel() =
BitRefLevel().set(skeleton_bit) | BitRefLevel().set(aggregate_bit);
simple->getBitRefLevelMask() = BitRefLevel().set();
simple->reSetUp(true);
CHKERR solveAdvection();
MoFEMFunctionReturn(0);
}
struct LevelSet::OpRhsDomain : public AssemblyDomainEleOp {
OpRhsDomain(const std::string field_name,
boost::shared_ptr<MatrixDouble> l_ptr,
boost::shared_ptr<MatrixDouble> l_dot_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr);
MoFEMErrorCode iNtegrate(EntData &data);
private:
boost::shared_ptr<MatrixDouble> lPtr;
boost::shared_ptr<MatrixDouble> lDotPtr;
boost::shared_ptr<MatrixDouble> velPtr;
};
struct LevelSet::OpLhsDomain : public AssemblyDomainEleOp {
OpLhsDomain(const std::string field_name,
boost::shared_ptr<MatrixDouble> vel_ptr);
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
private:
boost::shared_ptr<MatrixDouble> velPtr;
};
struct LevelSet::OpRhsSkeleton : public BoundaryEleOp {
OpRhsSkeleton(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr);
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<FaceSideEle>
sideFEPtr; ///< pointer to element to get data on edge/face sides
VectorDouble resSkelton;
};
struct LevelSet::OpLhsSkeleton : public BoundaryEleOp {
OpLhsSkeleton(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr);
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<FaceSideEle>
sideFEPtr; ///< pointer to element to get data on edge/face sides
MatrixDouble matSkeleton;
};
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
try {
// Create MoAB database
moab::Core moab_core;
moab::Interface &moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface &m_field = mofem_core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "LevelSet"));
LogManager::setLog("LevelSet");
MOFEM_LOG_TAG("LevelSet", "LevelSet");
LevelSet level_set(m_field);
CHKERR level_set.runProblem();
}
CATCH_ERRORS;
// finish work cleaning memory, getting statistics, etc.
MoFEM::Core::Finalize();
return 0;
}
MoFEMErrorCode LevelSet::readMesh() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
// get options from command line
CHKERR simple->getOptions();
// Only L2 field is set in this example. Two lines bellow forces simple
// interface to creat lower dimension (edge) elements, despite that fact that
// there is no field spanning on such elements. We need them for DG method.
simple->getAddSkeletonFE() = true;
simple->getAddBoundaryFE() = true;
// load mesh file
simple->getBitRefLevel() = BitRefLevel();
CHKERR simple->loadFile();
// Initialise bit ref levels
auto set_problem_bit = [&]() {
MoFEMFunctionBegin;
auto bit0 = BitRefLevel().set(start_bit);
BitRefLevel start_mask;
for (auto s = 0; s != start_bit; ++s)
start_mask[s] = true;
auto bit_mng = mField.getInterface<BitRefManager>();
Range level0;
CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(0),
BitRefLevel().set(), level0);
CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
// Set bits to build adjacencies between parents and children. That is used
// by simple interface.
simple->getBitAdjEnt() = BitRefLevel().set();
simple->getBitAdjParent() = BitRefLevel().set();
simple->getBitRefLevel() = BitRefLevel().set(current_bit);
simple->getBitRefLevelMask() = BitRefLevel().set();
#ifndef NDEBUG
if constexpr (debug) {
auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
CHKERR bit_mng->writeBitLevelByDim(
BitRefLevel().set(0), BitRefLevel().set(), FE_DIM,
(proc_str + "level_base.vtk").c_str(), "VTK", "");
}
#endif
MoFEMFunctionReturn(0);
};
CHKERR set_problem_bit();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode LevelSet::setUpProblem() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
// Scalar fields and vector field is tested. Add more fields, i.e. vector
// field if needed.
CHKERR simple->addDomainField("L", L2, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addDomainField("V", potential_velocity_space,
AINSWORTH_LEGENDRE_BASE, 1);
// set fields order, i.e. for most first cases order is sufficient.
CHKERR simple->setFieldOrder("L", 4);
CHKERR simple->setFieldOrder("V", 4);
// setup problem
CHKERR simple->setUp();
MoFEMFunctionReturn(0);
}
FTensor::Index<'i', SPACE_DIM> i;
FTensor::Index<'j', SPACE_DIM> j;
FTensor::Index<'k', SPACE_DIM> k;
LevelSet::OpRhsDomain::OpRhsDomain(const std::string field_name,
boost::shared_ptr<MatrixDouble> l_ptr,
boost::shared_ptr<MatrixDouble> l_dot_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr)
: AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
lPtr(l_ptr), lDotPtr(l_dot_ptr), velPtr(vel_ptr) {}
LevelSet::OpLhsDomain::OpLhsDomain(const std::string field_name,
boost::shared_ptr<MatrixDouble> vel_ptr)
: AssemblyDomainEleOp(field_name, field_name,
AssemblyDomainEleOp::OPROWCOL),
velPtr(vel_ptr) {
this->sYmm = false;
}
LevelSet::OpRhsSkeleton::OpRhsSkeleton(
boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
LevelSet::OpLhsSkeleton::OpLhsSkeleton(
boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr) {}
MoFEMErrorCode LevelSet::OpRhsDomain::iNtegrate(EntData &data) {
MoFEMFunctionBegin;
const auto nb_int_points = getGaussPts().size2();
const auto nb_dofs = data.getIndices().size();
const auto nb_base_func = data.getN().size2();
auto t_l = getFTensor2FromMat<DIM1, DIM2>(*lPtr);
auto t_l_dot = getFTensor2FromMat<DIM1, DIM2>(*lDotPtr);
auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velPtr);
auto t_base = data.getFTensor0N();
auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != nb_int_points; ++gg) {
const auto alpha = t_w * getMeasure();
FTensor::Tensor2<double, DIM1, DIM2> t_res0;
t_res0(I, J) = alpha * t_l_dot(I, J);
FTensor::Tensor3<double, SPACE_DIM, DIM1, DIM2> t_res1;
t_res1(i, I, J) = (alpha * t_l(I, J)) * t_vel(i);
++t_w;
++t_l;
++t_l_dot;
++t_vel;
auto &nf = this->locF;
auto t_nf = getFTensor2FromPtr<DIM1, DIM2>(&*nf.begin());
int rr = 0;
for (; rr != nb_dofs; ++rr) {
t_nf(I, J) += t_res0(I, J) * t_base - t_res1(i, I, J) * t_diff_base(i);
++t_base;
++t_diff_base;
++t_nf;
}
for (; rr < nb_base_func; ++rr) {
++t_base;
++t_diff_base;
}
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode LevelSet::OpLhsDomain::iNtegrate(EntData &row_data,
EntData &col_data) {
MoFEMFunctionBegin;
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_row_base = row_data.getFTensor0N();
auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
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) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_mat = getFTensor2FromPtr<DIM1, DIM2>(&mat(rr * DIM1, 0));
for (int cc = 0; cc != nb_col_dofs; ++cc) {
t_mat(I, J) +=
(beta * t_row_base - alpha * (t_row_diff_base(i) * t_vel(i))) *
t_col_base;
++t_col_base;
++t_mat;
}
++t_row_base;
++t_row_diff_base;
}
for (; rr < nb_base_func; ++rr) {
++t_row_base;
++t_row_diff_base;
}
++t_vel;
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode
LevelSet::OpRhsSkeleton::doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
// Collect data from side domain elements
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
auto not_side = [](auto s) {
return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
};
auto get_ntensor = [](auto &base_mat) {
return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
&*base_mat.data().begin());
};
if (in_the_loop > 0) {
// get normal of the face or edge
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
// gent number of DOFs on the right side.
const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
if (nb_rows) {
resSkelton.resize(nb_rows, false);
resSkelton.clear();
// get orientation of the local element edge
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)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Should be opposite sign");
#endif
// iterate the side cols
const auto nb_row_base_functions =
sideDataPtr->rowBaseSideMap[s0].size2();
auto t_w = getFTensor0IntegrationWeight();
auto arr_t_l = make_array(
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
auto arr_t_vel = make_array(
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())
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Inconsistent number of DOFs");
#endif
auto t_row_base = get_ntensor(sideDataPtr->rowBaseSideMap[s0]);
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
FTensor::Tensor1<double, SPACE_DIM> t_vel;
t_vel(i) = (arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
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];
FTensor::Tensor2<double, DIM1, DIM2> t_res;
t_res(I, J) = t_w * dot * l_upwind(I, J);
next();
++t_w;
auto t_res_skeleton =
getFTensor2FromPtr<DIM1, DIM2>(&*resSkelton.data().begin());
auto rr = 0;
for (; rr != nb_rows; ++rr) {
t_res_skeleton(I, J) += t_row_base * t_res(I, J);
++t_row_base;
++t_res_skeleton;
}
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
}
}
// assemble local operator vector to global vector
CHKERR ::VecSetValues(getTSf(),
sideDataPtr->indicesRowSideMap[s0].size(),
&*sideDataPtr->indicesRowSideMap[s0].begin(),
&*resSkelton.begin(), ADD_VALUES);
}
}
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode
LevelSet::OpLhsSkeleton::doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
// Collect data from side domain elements
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
auto not_side = [](auto s) {
return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
};
auto get_ntensor = [](auto &base_mat) {
return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
&*base_mat.data().begin());
};
if (in_the_loop > 0) {
// get normal of the face or edge
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
// gent number of DOFs on the right side.
const auto nb_rows = sideDataPtr->indicesRowSideMap[s0].size();
if (nb_rows) {
// get orientation of the local element edge
const auto opposite_s0 = not_side(s0);
const auto sense_row = sideDataPtr->senseMap[s0];
// iterate the side cols
const auto nb_row_base_functions =
sideDataPtr->rowBaseSideMap[s0].size2();
for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
// gent number of DOFs on the right side.
const auto nb_cols = sideDataPtr->indicesColSideMap[s1].size();
const auto sense_col = sideDataPtr->senseMap[s1];
// resize local element matrix
matSkeleton.resize(nb_rows, nb_cols, false);
matSkeleton.clear();
auto t_w = getFTensor0IntegrationWeight();
auto arr_t_vel = make_array(
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) {
FTensor::Tensor1<double, SPACE_DIM> t_vel;
t_vel(i) =
(arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
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];
FTensor::Tensor2<double, DIM1, DIM2> t_res;
t_res(I, J) = 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);
return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(ptr);
};
auto t_col_base =
get_ntensor(sideDataPtr->colBaseSideMap[s1], gg, 0);
auto t_mat_skeleton =
getFTensor2FromPtr<DIM1, DIM2>(&matSkeleton(rr * DIM1, 0));
FTensor::Tensor2<double, DIM1, DIM2> t_res_row;
t_res_row(I, J) = t_res(I, J) * t_row_base;
++t_row_base;
// iterate columns
for (size_t cc = 0; cc != nb_cols; ++cc) {
t_mat_skeleton(I, J) += t_res_row(I, J) * t_col_base;
++t_col_base;
++t_mat_skeleton;
}
}
}
for (; rr < nb_row_base_functions; ++rr) {
++t_row_base;
}
}
// assemble system
CHKERR ::MatSetValues(getTSB(),
sideDataPtr->indicesRowSideMap[s0].size(),
&*sideDataPtr->indicesRowSideMap[s0].begin(),
sideDataPtr->indicesColSideMap[s1].size(),
&*sideDataPtr->indicesColSideMap[s1].begin(),
&*matSkeleton.data().begin(), ADD_VALUES);
}
}
}
}
MoFEMFunctionReturn(0);
}
ForcesAndSourcesCore::UserDataOperator *
LevelSet::getZeroLevelVelOp(boost::shared_ptr<MatrixDouble> vel_ptr) {
auto get_parent_vel_this = [&]() {
auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
parent_fe_ptr->getOpPtrVector(), {potential_velocity_space});
parent_fe_ptr->getOpPtrVector().push_back(
new OpCalculateHcurlVectorCurl<potential_velocity_field_dim, SPACE_DIM>(
"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;
for (int l = 0; l <= nb_levels; ++l)
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
for (auto l = 1; l <= nb_levels; ++l) {
parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
BitRefLevel().set(0).flip(), this_fe_ptr,
BitRefLevel().set(0), BitRefLevel().set()));
}
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);
return new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
BitRefLevel().set(0).flip(), this_fe_ptr,
BitRefLevel().set(0), BitRefLevel().set());
}
MoFEMErrorCode LevelSet::pushOpDomain() {
MoFEMFunctionBegin;
auto pip = mField.getInterface<PipelineManager>(); // get interface to
// pipeline manager
pip->getOpDomainLhsPipeline().clear();
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(
current_bit);
};
pip->getDomainRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
current_bit);
};
auto l_ptr = boost::make_shared<MatrixDouble>();
auto l_dot_ptr = boost::make_shared<MatrixDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
pip->getOpDomainRhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
{L2});
pip->getOpDomainRhsPipeline().push_back(
new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
pip->getOpDomainRhsPipeline().push_back(
new OpCalculateTensor2FieldValuesDot<DIM1, DIM2>("L", l_dot_ptr));
pip->getOpDomainRhsPipeline().push_back(
new OpRhsDomain("L", l_ptr, l_dot_ptr, vel_ptr));
pip->getOpDomainLhsPipeline().push_back(getZeroLevelVelOp(vel_ptr));
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
{L2});
pip->getOpDomainLhsPipeline().push_back(new OpLhsDomain("L", vel_ptr));
MoFEMFunctionReturn(0);
}
boost::shared_ptr<FaceSideEle>
LevelSet::getSideFE(boost::shared_ptr<SideData> side_data_ptr) {
auto simple = mField.getInterface<Simple>();
auto l_ptr = boost::make_shared<MatrixDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
struct OpSideData : public FaceSideEleOp {
OpSideData(boost::shared_ptr<SideData> side_data_ptr)
: FaceSideEleOp("L", "L", FaceSideEleOp::OPROWCOL),
sideDataPtr(side_data_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
doEntities[t] = true;
sYmm = false;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
MoFEMFunctionBegin;
if ((CN::Dimension(row_type) == FE_DIM) &&
(CN::Dimension(col_type) == FE_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 {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
}
MoFEMFunctionReturn(0);
};
private:
boost::shared_ptr<SideData> sideDataPtr;
};
struct OpSideDataOnParent : public DomainEleOp {
OpSideDataOnParent(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<MatrixDouble> l_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr)
: DomainEleOp("L", "L", DomainEleOp::OPROWCOL),
sideDataPtr(side_data_ptr), lPtr(l_ptr), velPtr(vel_ptr) {
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
for (auto t = moab::CN::TypeDimensionMap[FE_DIM].first;
t <= moab::CN::TypeDimensionMap[FE_DIM].second; ++t)
doEntities[t] = true;
sYmm = false;
}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type, EntData &row_data,
EntData &col_data) {
MoFEMFunctionBegin;
if ((CN::Dimension(row_type) == FE_DIM) &&
(CN::Dimension(col_type) == FE_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].size2() !=
(sideDataPtr->velMat)[nb_in_loop].size2())
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of integaration pts %d != %d",
(sideDataPtr->lVec)[nb_in_loop].size2(),
(sideDataPtr->velMat)[nb_in_loop].size2());
if ((sideDataPtr->velMat)[nb_in_loop].size1() != SPACE_DIM)
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"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()) {
SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of integration pt %d != %d",
sideDataPtr->rowBaseSideMap[0].size1(),
sideDataPtr->rowBaseSideMap[1].size1());
}
if (sideDataPtr->colBaseSideMap[0].size1() !=
sideDataPtr->colBaseSideMap[1].size1()) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Wrong number of integration pt");
}
#endif
}
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
}
MoFEMFunctionReturn(0);
};
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<MatrixDouble> lPtr;
boost::shared_ptr<MatrixDouble> velPtr;
};
// Calculate fields on param mesh bit element
auto get_parent_this = [&]() {
auto parent_fe_ptr = boost::make_shared<DomianParentEle>(mField);
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
parent_fe_ptr->getOpPtrVector(), {L2});
parent_fe_ptr->getOpPtrVector().push_back(
new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
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;
for (int l = 0; l <= nb_levels; ++l)
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
for (auto l = 1; l <= nb_levels; ++l) {
parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
BitRefLevel().set(current_bit).flip(), this_fe_ptr,
BitRefLevel().set(current_bit), BitRefLevel().set()));
}
return parents_elems_ptr_vec[0];
};
// Create aliased shared pointers, all elements are destroyed if side_fe_ptr
// is destroyed
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(getZeroLevelVelOp(vel_ptr));
side_fe_ptr->getOpPtrVector().push_back(
new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
BitRefLevel().set(current_bit).flip(), this_fe_ptr,
BitRefLevel().set(current_bit), BitRefLevel().set()));
return side_fe_ptr;
};
return get_side_fe_ptr();
};
MoFEMErrorCode LevelSet::pushOpSkeleton() {
MoFEMFunctionBegin;
auto pip = mField.getInterface<PipelineManager>(); // get interface to
pip->getOpSkeletonLhsPipeline().clear();
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(
skeleton_bit);
};
pip->getSkeletonRhsFE()->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
skeleton_bit);
};
auto side_data_ptr = boost::make_shared<SideData>();
auto side_fe_ptr = getSideFE(side_data_ptr);
pip->getOpSkeletonRhsPipeline().push_back(
new OpRhsSkeleton(side_data_ptr, side_fe_ptr));
pip->getOpSkeletonLhsPipeline().push_back(
new OpLhsSkeleton(side_data_ptr, side_fe_ptr));
MoFEMFunctionReturn(0);
}
std::tuple<double, Tag> LevelSet::evaluateError() {
struct OpErrorSkel : BoundaryEleOp {
OpErrorSkel(boost::shared_ptr<FaceSideEle> side_fe_ptr,
boost::shared_ptr<SideData> side_data_ptr,
SmartPetscObj<Vec> error_sum_ptr, Tag th_error)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
sideFEPtr(side_fe_ptr), sideDataPtr(side_data_ptr),
errorSumPtr(error_sum_ptr), thError(th_error) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
// Collect data from side domain elements
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
auto not_side = [](auto s) {
return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
};
auto nb_gauss_pts = getGaussPts().size2();
for (auto s : {LEFT_SIDE, RIGHT_SIDE}) {
auto arr_t_l = make_array(
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
auto arr_t_vel = make_array(
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) {
FTensor::Tensor2<double, DIM1, DIM2> t_diff;
t_diff(I, J) = arr_t_l[LEFT_SIDE](I, J) - arr_t_l[RIGHT_SIDE](I, J);
e += t_w * getMeasure() * t_diff(I, J) * t_diff(I, J);
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);
};
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<FaceSideEle> sideFEPtr;
boost::shared_ptr<SideData> sideDataPtr;
SmartPetscObj<Vec> errorSumPtr;
Tag thError;
};
auto simple = mField.getInterface<Simple>();
auto error_sum_ptr = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
Tag th_error;
double def_val = 0;
CHKERR mField.get_moab().tag_get_handle("Error", 1, MB_TYPE_DOUBLE, th_error,
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_val);
auto clear_tags = [&]() {
MoFEMFunctionBegin;
Range fe_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, FE_DIM, fe_ents);
double zero;
CHKERR mField.get_moab().tag_clear_data(th_error, fe_ents, &zero);
MoFEMFunctionReturn(0);
};
auto evaluate_error = [&]() {
MoFEMFunctionBegin;
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>();
auto side_fe_ptr = getSideFE(side_data_ptr);
skel_fe->getOpPtrVector().push_back(
new OpErrorSkel(side_fe_ptr, side_data_ptr, error_sum_ptr, th_error));
auto simple = mField.getInterface<Simple>();
skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
skeleton_bit);
};
CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
simple->getSkeletonFEName(), skel_fe);
MoFEMFunctionReturn(0);
};
auto assemble_and_sum = [](auto vec) {
CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
double sum;
CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
return sum;
};
auto propagate_error_to_parents = [&]() {
MoFEMFunctionBegin;
auto &moab = mField.get_moab();
auto fe_ptr = boost::make_shared<FEMethod>();
fe_ptr->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
current_bit);
};
fe_ptr->preProcessHook = []() { return 0; };
fe_ptr->postProcessHook = []() { return 0; };
fe_ptr->operatorHook = [&]() {
MoFEMFunctionBegin;
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);
boost::function<MoFEMErrorCode(EntityHandle, double)> add_error =
[&](auto fe_ent, auto error) {
MoFEMFunctionBegin;
double *e_ptr;
CHKERR moab.tag_get_by_ptr(th_error, &fe_ent, 1,
(const void **)&e_ptr);
(*e_ptr) += error;
EntityHandle parent;
CHKERR moab.tag_get_data(th_parent, &fe_ent, 1, &parent);
if (parent != fe_ent && parent)
CHKERR add_error(parent, *e_ptr);
MoFEMFunctionReturn(0);
};
CHKERR add_error(parent, error);
MoFEMFunctionReturn(0);
};
CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
fe_ptr);
MoFEMFunctionReturn(0);
};
CHK_THROW_MESSAGE(clear_tags(), "clear error tags");
CHK_THROW_MESSAGE(evaluate_error(), "evaluate error");
CHK_THROW_MESSAGE(propagate_error_to_parents(), "propagate error");
return std::make_tuple(assemble_and_sum(error_sum_ptr), th_error);
}
/**
* @brief test side element
*
* Check consistency between volume and skeleton integral
*
* @return MoFEMErrorCode
*/
MoFEMErrorCode LevelSet::testSideFE() {
MoFEMFunctionBegin;
/**
* @brief calculate volume
*
*/
struct DivergenceVol : public DomainEleOp {
DivergenceVol(boost::shared_ptr<MatrixDouble> l_ptr,
boost::shared_ptr<MatrixDouble> vel_ptr,
SmartPetscObj<Vec> div_vec)
: DomainEleOp("L", DomainEleOp::OPROW), lPtr(l_ptr), velPtr(vel_ptr),
divVec(div_vec) {}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
const auto nb_dofs = data.getIndices().size();
if (nb_dofs) {
const auto nb_gauss_pts = getGaussPts().size2();
const auto t_w = getFTensor0IntegrationWeight();
auto t_diff = data.getFTensor1DiffN<SPACE_DIM>();
auto t_l = getFTensor2FromMat<DIM1, DIM2>(*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(I, I) * (t_diff(i) * t_vel(i));
++t_diff;
}
++t_w;
++t_l;
++t_vel;
}
CHKERR VecSetValue(divVec, 0, div, ADD_VALUES);
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<MatrixDouble> lPtr;
boost::shared_ptr<MatrixDouble> velPtr;
SmartPetscObj<Vec> divVec;
};
/**
* @brief calculate skeleton integral
*
*/
struct DivergenceSkeleton : public BoundaryEleOp {
DivergenceSkeleton(boost::shared_ptr<SideData> side_data_ptr,
boost::shared_ptr<FaceSideEle> side_fe_ptr,
SmartPetscObj<Vec> div_vec)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
sideDataPtr(side_data_ptr), sideFEPtr(side_fe_ptr), divVec(div_vec) {}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
auto get_ntensor = [](auto &base_mat) {
return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
&*base_mat.data().begin());
};
auto not_side = [](auto s) {
return s == LEFT_SIDE ? RIGHT_SIDE : LEFT_SIDE;
};
// Collect data from side domain elements
CHKERR loopSideFaces("dFE", *sideFEPtr);
const auto in_the_loop =
sideFEPtr->nInTheLoop; // return number of elements on the side
auto t_normal = getFTensor1Normal();
const auto nb_gauss_pts = getGaussPts().size2();
for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
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 = make_array(
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[LEFT_SIDE]),
getFTensor2FromMat<DIM1, DIM2>(sideDataPtr->lVec[RIGHT_SIDE]));
auto arr_t_vel = make_array(
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) {
FTensor::Tensor1<double, SPACE_DIM> t_vel;
t_vel(i) =
(arr_t_vel[LEFT_SIDE](i) + arr_t_vel[RIGHT_SIDE](i)) / 2.;
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]; //< assume that field is continues,
// initialisation field has to be smooth
// and exactly approximated by approx
// base
auto res = t_w * l_upwind(I, I) * 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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<SideData> sideDataPtr;
boost::shared_ptr<FaceSideEle> sideFEPtr;
boost::shared_ptr<MatrixDouble> velPtr;
SmartPetscObj<Vec> divVec;
};
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 div_vol_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
auto div_skel_vec = createVectorMPI(mField.get_comm(), PETSC_DECIDE, 1);
auto l_ptr = boost::make_shared<MatrixDouble>();
auto vel_ptr = boost::make_shared<MatrixDouble>();
auto side_data_ptr = boost::make_shared<SideData>();
auto side_fe_ptr = getSideFE(side_data_ptr);
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(vol_fe->getOpPtrVector(),
{L2});
vol_fe->getOpPtrVector().push_back(
new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_ptr));
vol_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
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));
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
/**
* Set up problem such that gradient of level set field is orthogonal to
* velocity field. Then volume and skeleton integral should yield the same
* value.
*/
CHKERR initialiseFieldVelocity(
[](double x, double y, double) { return x - y; });
CHKERR initialiseFieldLevelSet(
[](double x, double y, double) { return x - y; });
vol_fe->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
current_bit);
};
skel_fe->exeTestHook = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
skeleton_bit);
};
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_fe);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getSkeletonFEName(), skel_fe);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), skel_fe);
auto assemble_and_sum = [](auto vec) {
CHK_THROW_MESSAGE(VecAssemblyBegin(vec), "assemble");
CHK_THROW_MESSAGE(VecAssemblyEnd(vec), "assemble");
double sum;
CHK_THROW_MESSAGE(VecSum(vec, &sum), "assemble");
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;
if (eps > eps_err)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"No consistency between skeleton integral and volume integral");
MoFEMFunctionReturn(0);
};
MoFEMErrorCode LevelSet::testOp() {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto opt = mField.getInterface<OperatorsTester>(); // get interface to
// OperatorsTester
auto pip = mField.getInterface<PipelineManager>(); // get interface to
// pipeline manager
CHKERR pushOpDomain();
CHKERR pushOpSkeleton();
auto post_proc = [&](auto dm, auto f_res, auto out_name) {
MoFEMFunctionBegin;
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
if constexpr (DIM1 == 1 && DIM2 == 1) {
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
auto l_vec = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", l_vec, f_res));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{},
{}, {})
);
}
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
SmartPetscObj<Vec>(), diff_x, 0, 1, eps);
if constexpr (debug) {
// Example how to plot direction in direction diff_x. If instead
// directionalCentralFiniteDifference(...) diff_res is used, then error
// on directive is plotted.
CHKERR post_proc(simple->getDM(), diff_res, "tangent_op_error.h5m");
}
// Calculate norm of difference between directive calculated from finite
// difference, and tangent matrix.
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C("LevelSet", Sev::inform,
"Test consistency of tangent matrix %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err)
SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Norm of directional derivative too large err = %3.4e", fnorm);
MoFEMFunctionReturn(0);
};
CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
pip->getDomainRhsFE());
CHKERR test_domain_ops(simple->getSkeletonFEName(), pip->getSkeletonLhsFE(),
pip->getSkeletonRhsFE());
MoFEMFunctionReturn(0);
};
MoFEMErrorCode LevelSet::initialiseFieldLevelSet(
boost::function<double(double, double, double)> level_fun) {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>(); // get interface to
// pipeline manager
auto prb_mng = mField.getInterface<ProblemsManager>();
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; });
auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "LEVELSET_POJECTION");
CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
CHKERR DMSetUp(sub_dm);
BitRefLevel remove_mask = BitRefLevel().set(current_bit);
remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
CHKERR prb_mng->removeDofsOnEntities("LEVELSET_POJECTION", "L",
BitRefLevel().set(), remove_mask);
auto test_bit_ele = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
current_bit);
};
pip->getDomainLhsFE()->exeTestHook = test_bit_ele;
pip->getDomainRhsFE()->exeTestHook = test_bit_ele;
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
{L2});
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
{L2});
pip->getOpDomainLhsPipeline().push_back(new OpMassLL("L", "L"));
pip->getOpDomainRhsPipeline().push_back(new OpSourceL("L", level_fun));
CHKERR mField.getInterface<FieldBlas>()->setField(0, "L");
auto ksp = pip->createKSP(sub_dm);
CHKERR KSPSetDM(ksp, sub_dm);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSetUp(ksp);
auto L = createDMVector(sub_dm);
auto F = vectorDuplicate(L);
CHKERR KSPSolve(ksp, F, L);
CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
auto [error, th_error] = evaluateError();
MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
#ifndef NDEBUG
auto fe_meshset =
mField.get_finite_element_meshset(simple->getDomainFEName());
std::vector<Tag> tags{th_error};
CHKERR mField.get_moab().write_file("error.h5m", "MOAB",
"PARALLEL=WRITE_PART", &fe_meshset, 1,
&*tags.begin(), tags.size());
#endif
auto post_proc = [&](auto dm, auto out_name, auto th_error) {
MoFEMFunctionBegin;
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->setTagsToTransfer({th_error});
post_proc_fe->exeTestHook = test_bit_ele;
if constexpr (DIM1 == 1 && DIM2 == 1) {
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
auto l_vec = boost::make_shared<VectorDouble>();
auto l_grad_mat = boost::make_shared<MatrixDouble>();
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {L2});
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", l_vec));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{{"GradL", l_grad_mat}},
{}, {})
);
}
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
MoFEMFunctionReturn(0);
};
if constexpr (debug)
CHKERR post_proc(sub_dm, "initial_level_set.h5m", th_error);
swap_fe();
MoFEMFunctionReturn(0);
}
MoFEMErrorCode LevelSet::initialiseFieldVelocity(
boost::function<double(double, double, double)> vel_fun) {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>(); // get interface to
// pipeline manager
auto prb_mng = mField.getInterface<ProblemsManager>();
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; });
auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "VELOCITY_PROJECTION");
CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
CHKERR DMMoFEMAddSubFieldRow(sub_dm, "V");
CHKERR DMSetUp(sub_dm);
// Velocities are calculated only on corse mesh
BitRefLevel remove_mask = BitRefLevel().set(0);
remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
CHKERR prb_mng->removeDofsOnEntities("VELOCITY_PROJECTION", "V",
BitRefLevel().set(), remove_mask);
auto test_bit = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(0);
};
pip->getDomainLhsFE()->exeTestHook = test_bit;
pip->getDomainRhsFE()->exeTestHook = test_bit;
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainLhsPipeline(),
{potential_velocity_space});
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(pip->getOpDomainRhsPipeline(),
{potential_velocity_space});
pip->getOpDomainLhsPipeline().push_back(new OpMassVV("V", "V"));
pip->getOpDomainRhsPipeline().push_back(new OpSourceV("V", vel_fun));
auto ksp = pip->createKSP(sub_dm);
CHKERR KSPSetDM(ksp, sub_dm);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSetUp(ksp);
auto L = createDMVector(sub_dm);
auto F = vectorDuplicate(L);
CHKERR KSPSolve(ksp, F, L);
CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
auto post_proc = [&](auto dm, auto out_name) {
MoFEMFunctionBegin;
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->exeTestHook = test_bit;
if constexpr (FE_DIM == 2) {
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {potential_velocity_space});
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
auto potential_vec = boost::make_shared<VectorDouble>();
auto velocity_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("V", potential_vec));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHcurlVectorCurl<potential_velocity_field_dim,
SPACE_DIM>("V", velocity_mat));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"VelocityPotential", potential_vec}},
{{"Velocity", velocity_mat}},
{}, {})
);
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"3d case not implemented");
}
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
MoFEMFunctionReturn(0);
};
if constexpr (debug)
CHKERR post_proc(sub_dm, "initial_velocity_potential.h5m");
swap_fe();
MoFEMFunctionReturn(0);
}
LevelSet *level_set_raw_ptr = nullptr;
MoFEM::TsCtx *ts_ctx;
MoFEMErrorCode LevelSet::solveAdvection() {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>(); // get interface to
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR pushOpDomain();
CHKERR pushOpSkeleton();
auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "ADVECTION");
CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
CHKERR DMMoFEMAddElement(sub_dm, simple->getSkeletonFEName());
CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
CHKERR DMSetUp(sub_dm);
BitRefLevel remove_mask = BitRefLevel().set(current_bit);
remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
CHKERR prb_mng->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
remove_mask);
auto add_post_proc_fe = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
Tag th_error;
double def_val = 0;
CHKERR mField.get_moab().tag_get_handle(
"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(
current_bit);
};
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
auto vel_ptr = boost::make_shared<MatrixDouble>();
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {L2});
post_proc_fe->getOpPtrVector().push_back(getZeroLevelVelOp(vel_ptr));
if constexpr (DIM1 == 1 && DIM2 == 1) {
auto l_vec = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", l_vec));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
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 = [&]() {
MoFEMFunctionBegin;
if (!post_proc_fe)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Null pointer for post proc element");
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
CHKERR post_proc_fe->writeFile(
"level_set_" +
boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
MoFEMFunctionReturn(0);
};
boost::shared_ptr<FEMethod> null;
DMMoFEMTSSetMonitor(sub_dm, ts, simple->getDomainFEName(), monitor_ptr,
null, null);
return monitor_ptr;
};
auto ts = pip->createTSIM(sub_dm);
auto set_solution = [&](auto ts) {
MoFEMFunctionBegin;
auto D = createDMVector(sub_dm);
CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR TSSetSolution(ts, D);
MoFEMFunctionReturn(0);
};
CHKERR set_solution(ts);
auto monitor_pt = set_time_monitor(sub_dm, ts);
CHKERR TSSetFromOptions(ts);
auto B = createDMMatrix(sub_dm);
CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
level_set_raw_ptr = this;
CHKERR TSSetUp(ts);
auto ts_pre_step = [](TS ts) {
auto &m_field = level_set_raw_ptr->mField;
auto simple = m_field.getInterface<Simple>();
auto bit_mng = m_field.getInterface<BitRefManager>();
MoFEMFunctionBegin;
auto [error, th_error] = level_set_raw_ptr->evaluateError();
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) {
MoFEMFunctionBegin;
DM dm;
CHKERR TSGetDM(ts, &dm);
auto prb_ptr = getProblemPtr(dm);
auto x = createDMVector(dm);
CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG("LevelSet", Sev::inform)
<< "Problem " << prb_ptr->getName() << " solution vector norm "
<< get_norm(x);
CHKERR TSSetSolution(ts, x);
MoFEMFunctionReturn(0);
};
auto refine_and_project = [&](auto ts) {
MoFEMFunctionBegin;
CHKERR level_set_raw_ptr->refineMesh(
WrapperClassErrorProjection(level_set_raw_ptr->maxPtr));
simple->getBitRefLevel() = BitRefLevel().set(skeleton_bit) |
BitRefLevel().set(aggregate_bit) |
BitRefLevel().set(aggregate_projection_bit);
simple->reSetUp(true);
DM dm;
CHKERR TSGetDM(ts, &dm);
CHKERR DMSubDMSetUp_MoFEM(dm);
BitRefLevel remove_mask = BitRefLevel().set(current_bit);
remove_mask
.flip(); // DOFs which are not on bit_domain_ele should be removed
CHKERR level_set_raw_ptr->mField.getInterface<ProblemsManager>()
->removeDofsOnEntities("ADVECTION", "L", BitRefLevel().set(),
remove_mask);
MoFEMFunctionReturn(0);
};
auto ts_reset_theta = [&](auto ts) {
MoFEMFunctionBegin;
DM dm;
CHKERR TSGetDM(ts, &dm);
// FIXME: Look into vec-5 how to transfer internal theta method variables
CHKERR TSReset(ts);
CHKERR TSSetUp(ts);
CHKERR level_set_raw_ptr->dgProjection(projection_bit);
CHKERR set_solution(ts);
auto B = createDMMatrix(dm);
CHKERR TSSetIJacobian(ts, B, B, TsSetIJacobian, nullptr);
MoFEMFunctionReturn(0);
};
CHKERR refine_and_project(ts);
CHKERR ts_reset_theta(ts);
MoFEMFunctionReturn(0);
};
auto ts_post_step = [](TS ts) { return 0; };
CHKERR TSSetPreStep(ts, ts_pre_step);
CHKERR TSSetPostStep(ts, ts_post_step);
CHKERR TSSolve(ts, NULL);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode LevelSet::refineMesh(WrapperClass &&wp) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto bit_mng = mField.getInterface<BitRefManager>();
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
auto set_bit = [](auto l) { return BitRefLevel().set(l); };
auto save_range = [&](const std::string name, const Range &r) {
MoFEMFunctionBegin;
auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
CHKERR mField.get_moab().add_entities(*meshset_ptr, r);
CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
meshset_ptr->get_ptr(), 1);
MoFEMFunctionReturn(0);
};
// select domain elements to refine by threshold
auto get_refined_elements_meshset = [&](auto bit, auto mask) {
Range fe_ents;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit, mask, FE_DIM, fe_ents);
Tag th_error;
CHK_MOAB_THROW(mField.get_moab().tag_get_handle("Error", th_error),
"get error handle");
std::vector<double> errors(fe_ents.size());
CHK_MOAB_THROW(
mField.get_moab().tag_get_data(th_error, fe_ents, &*errors.begin()),
"get tag data");
auto it = std::max_element(errors.begin(), errors.end());
double max;
MPI_Allreduce(&*it, &max, 1, MPI_DOUBLE, MPI_MAX, mField.get_comm());
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;
}
Range ents;
ents.insert_list(fe_to_refine.begin(), fe_to_refine.end());
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
ents, nullptr, NOISY);
auto get_neighbours_by_bridge_vertices = [&](auto &&ents) {
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts, true);
CHKERR mField.get_moab().get_adjacencies(verts, FE_DIM, false, ents,
moab::Interface::UNION);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
return ents;
};
ents = get_neighbours_by_bridge_vertices(ents);
#ifndef NDEBUG
if (debug) {
auto meshset_ptr = get_temp_meshset_ptr(mField.get_moab());
CHK_MOAB_THROW(mField.get_moab().add_entities(*meshset_ptr, ents),
"add entities to meshset");
CHKERR mField.get_moab().write_file(
(proc_str + "_fe_to_refine.vtk").c_str(), "VTK", "",
meshset_ptr->get_ptr(), 1);
}
#endif
return ents;
};
// refine elements, and set bit ref level
auto refine_mesh = [&](auto l, auto &&fe_to_refine) {
Skinner skin(&mField.get_moab());
MoFEMFunctionBegin;
// get entities in "l-1" level
Range level_ents;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
set_bit(start_bit + l - 1), BitRefLevel().set(), FE_DIM, level_ents);
// select entities to refine
fe_to_refine = intersect(level_ents, fe_to_refine);
// select entities not to refine
level_ents = subtract(level_ents, fe_to_refine);
// for entities to refine get children, i.e. redlined entities
Range fe_to_refine_children;
bit_mng->updateRangeByChildren(fe_to_refine, fe_to_refine_children);
// add entities to to level "l"
fe_to_refine_children = fe_to_refine_children.subset_by_dimension(FE_DIM);
level_ents.merge(fe_to_refine_children);
auto fix_neighbour_level = [&](auto ll) {
MoFEMFunctionBegin;
// filter entities on level ll
auto level_ll = level_ents;
CHKERR bit_mng->filterEntitiesByRefLevel(set_bit(ll), BitRefLevel().set(),
level_ll);
// find skin of ll level
Range skin_edges;
CHKERR skin.find_skin(0, level_ll, false, skin_edges);
// get parents of skin of level ll
Range skin_parents;
for (auto lll = 0; lll <= ll; ++lll) {
CHKERR bit_mng->updateRangeByParent(skin_edges, skin_parents);
}
// filter parents on level ll - 1
BitRefLevel bad_bit;
for (auto lll = 0; lll <= ll - 2; ++lll) {
bad_bit[lll] = true;
}
// get adjacents to parents
Range skin_adj_ents;
CHKERR mField.get_moab().get_adjacencies(
skin_parents, FE_DIM, false, skin_adj_ents, moab::Interface::UNION);
CHKERR bit_mng->filterEntitiesByRefLevel(bad_bit, BitRefLevel().set(),
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);
}
MoFEMFunctionReturn(0);
};
CHKERR fix_neighbour_level(l);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
level_ents);
// get lower dimension entities for level "l"
for (auto d = 0; d != FE_DIM; ++d) {
if (d == 0) {
CHKERR mField.get_moab().get_connectivity(
level_ents.subset_by_dimension(FE_DIM), level_ents, true);
} else {
CHKERR mField.get_moab().get_adjacencies(
level_ents.subset_by_dimension(FE_DIM), d, false, level_ents,
moab::Interface::UNION);
}
}
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
level_ents);
// set bit ref level to level entities
CHKERR bit_mng->setNthBitRefLevel(start_bit + l, false);
CHKERR bit_mng->setNthBitRefLevel(level_ents, start_bit + l, true);
#ifndef NDEBUG
auto proc_str = boost::lexical_cast<std::string>(mField.get_comm_rank());
CHKERR bit_mng->writeBitLevelByDim(
set_bit(start_bit + l), BitRefLevel().set(), FE_DIM,
(boost::lexical_cast<std::string>(l) + "_" + proc_str + "_ref_mesh.vtk")
.c_str(),
"VTK", "");
#endif
MoFEMFunctionReturn(0);
};
// set skeleton
auto set_skelton_bit = [&](auto l) {
MoFEMFunctionBegin;
// get entities of dim-1 on level "l"
Range level_edges;
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
set_bit(start_bit + l), BitRefLevel().set(), FE_DIM - 1, level_edges);
// get parent of entities of level "l"
Range level_edges_parents;
CHKERR bit_mng->updateRangeByParent(level_edges, level_edges_parents);
level_edges_parents = level_edges_parents.subset_by_dimension(FE_DIM - 1);
CHKERR bit_mng->filterEntitiesByRefLevel(
set_bit(start_bit + l), BitRefLevel().set(), level_edges_parents);
// skeleton entities which do not have parents
auto parent_skeleton = intersect(level_edges, level_edges_parents);
auto skeleton = subtract(level_edges, level_edges_parents);
// add adjacent domain entities
CHKERR mField.get_moab().get_adjacencies(unite(parent_skeleton, skeleton),
FE_DIM, false, skeleton,
moab::Interface::UNION);
// set levels
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(skeleton);
CHKERR bit_mng->setNthBitRefLevel(skeleton_bit, false);
CHKERR bit_mng->setNthBitRefLevel(skeleton, skeleton_bit, true);
#ifndef NDEBUG
CHKERR bit_mng->writeBitLevel(
set_bit(skeleton_bit), BitRefLevel().set(),
(boost::lexical_cast<std::string>(l) + "_" + proc_str + "_skeleton.vtk")
.c_str(),
"VTK", "");
#endif
MoFEMFunctionReturn(0);
};
// Reset bit sand set old current and aggregate bits as projection bits
Range level0_current;
CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(current_bit),
BitRefLevel().set(), level0_current);
Range level0_aggregate;
CHKERR bit_mng->getEntitiesByRefLevel(BitRefLevel().set(aggregate_bit),
BitRefLevel().set(), level0_aggregate);
BitRefLevel start_mask;
for (auto s = 0; s != start_bit; ++s)
start_mask[s] = true;
CHKERR bit_mng->lambdaBitRefLevel(
[&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
CHKERR bit_mng->setNthBitRefLevel(level0_current, projection_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0_aggregate, aggregate_projection_bit,
true);
// Set zero bit ref level
Range level0;
CHKERR bit_mng->getEntitiesByRefLevel(set_bit(0), BitRefLevel().set(),
level0);
CHKERR bit_mng->setNthBitRefLevel(level0, start_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0, current_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0, aggregate_bit, true);
CHKERR bit_mng->setNthBitRefLevel(level0, skeleton_bit, true);
CHKERR wp.setBits(*this, 0);
CHKERR wp.runCalcs(*this, 0);
for (auto l = 0; l != nb_levels; ++l) {
MOFEM_LOG("WORLD", Sev::inform) << "Process level: " << l;
CHKERR refine_mesh(l + 1, get_refined_elements_meshset(
set_bit(start_bit + l), BitRefLevel().set()));
CHKERR set_skelton_bit(l + 1);
CHKERR wp.setAggregateBit(*this, l + 1);
CHKERR wp.setBits(*this, l + 1);
CHKERR wp.runCalcs(*this, l + 1);
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode LevelSet::dgProjection(const int projection_bit) {
MoFEMFunctionBegin;
// get operators tester
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>(); // get interface to
auto bit_mng = mField.getInterface<BitRefManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
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; };
auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "DG_PROJECTION");
CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
CHKERR DMMoFEMAddSubFieldRow(sub_dm, "L");
CHKERR DMSetUp(sub_dm);
Range current_ents; // ents used to do calculations
CHKERR bit_mng->getEntitiesByDimAndRefLevel(BitRefLevel().set(current_bit),
BitRefLevel().set(), FE_DIM,
current_ents);
Range prj_ents; // ents from which data are projected
CHKERR bit_mng->getEntitiesByDimAndRefLevel(
BitRefLevel().set(projection_bit), BitRefLevel().set(), FE_DIM, prj_ents);
for (auto l = 0; l != nb_levels; ++l) {
CHKERR bit_mng->updateRangeByParent(prj_ents, prj_ents);
}
current_ents = subtract(
current_ents, prj_ents); // only restric to entities needed projection
auto test_mesh_bit = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
current_bit);
};
auto test_prj_bit = [&](FEMethod *fe_ptr) {
return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
projection_bit);
};
auto test_current_bit = [&](FEMethod *fe_ptr) {
return current_ents.find(fe_ptr->getFEEntityHandle()) != current_ents.end();
};
lhs_fe->exeTestHook =
test_mesh_bit; // that element only is run when current bit is set
rhs_fe_prj->exeTestHook =
test_prj_bit; // that element is run only when projection bit is set
rhs_fe_current->exeTestHook =
test_current_bit; // that element is only run when current bit is set
BitRefLevel remove_mask = BitRefLevel().set(current_bit);
remove_mask.flip(); // DOFs which are not on bit_domain_ele should be removed
CHKERR prb_mng->removeDofsOnEntities(
"DG_PROJECTION", "L", BitRefLevel().set(), remove_mask, nullptr, 0,
MAX_DOFS_ON_ENTITY, 0, 100, NOISY,
true); // remove all DOFs which are not
// on current bit. This case works for L2 space
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(lhs_fe->getOpPtrVector(),
{L2});
lhs_fe->getOpPtrVector().push_back(
new OpMassLL("L", "L")); // Assemble projection matrix
// This assumes that projection mesh is finer, current mesh is coarsened.
auto set_prj_from_child = [&](auto rhs_fe_prj) {
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
rhs_fe_prj->getOpPtrVector(), {L2});
// Evaluate field value on projection mesh
auto l_vec = boost::make_shared<MatrixDouble>();
rhs_fe_prj->getOpPtrVector().push_back(
new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_vec));
// This element is used to assemble
auto get_parent_this = [&]() {
auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
fe_parent_this->getOpPtrVector().push_back(
new OpScalarFieldL("L", l_vec));
return fe_parent_this;
};
// Create levels of parent elements, until current element is reached, and
// then assemble.
auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
for (int l = 0; l <= nb_levels; ++l)
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
for (auto l = 1; l <= nb_levels; ++l) {
parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
BitRefLevel().set(current_bit).flip(), this_fe_ptr,
BitRefLevel().set(current_bit),
BitRefLevel().set()));
}
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(
new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
BitRefLevel().set(current_bit).flip(), this_fe_ptr,
BitRefLevel().set(current_bit), BitRefLevel().set()));
};
// This assumed that current mesh is refined, and projection mesh is coarser
auto set_prj_from_parent = [&](auto rhs_fe_current) {
// Evaluate field value on projection mesh
auto l_vec = boost::make_shared<MatrixDouble>();
// Evaluate field on coarser element
auto get_parent_this = [&]() {
auto fe_parent_this = boost::make_shared<DomianParentEle>(mField);
fe_parent_this->getOpPtrVector().push_back(
new OpCalculateTensor2FieldValues<DIM1, DIM2>("L", l_vec));
return fe_parent_this;
};
// Create stack of evaluation on parent elements
auto get_parents_fe_ptr = [&](auto this_fe_ptr) {
std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
for (int l = 0; l <= nb_levels; ++l)
parents_elems_ptr_vec.emplace_back(
boost::make_shared<DomianParentEle>(mField));
for (auto l = 1; l <= nb_levels; ++l) {
parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
new OpRunParent(parents_elems_ptr_vec[l], BitRefLevel().set(),
BitRefLevel().set(projection_bit).flip(),
this_fe_ptr, BitRefLevel().set(projection_bit),
BitRefLevel().set()));
}
return parents_elems_ptr_vec[0];
};
auto this_fe_ptr = get_parent_this();
auto parent_fe_ptr = get_parents_fe_ptr(this_fe_ptr);
auto reset_op_ptr = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
reset_op_ptr->doWorkRhsHook = [&](DataOperator *op_ptr, int, EntityType,
EntData &) {
l_vec->resize(static_cast<DomainEleOp *>(op_ptr)->getGaussPts().size2(),
false);
l_vec->clear();
return 0;
};
rhs_fe_current->getOpPtrVector().push_back(reset_op_ptr);
rhs_fe_current->getOpPtrVector().push_back(
new OpRunParent(parent_fe_ptr, BitRefLevel().set(),
BitRefLevel().set(projection_bit).flip(), this_fe_ptr,
BitRefLevel(), BitRefLevel()));
// At the end assemble of current finite element
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;
getDMKspCtx(sub_dm)->clearLoops();
CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(),
lhs_fe, null_fe, null_fe);
CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), rhs_fe_prj,
null_fe, null_fe);
CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(),
rhs_fe_current, null_fe, null_fe);
auto ksp = MoFEM::createKSP(mField.get_comm());
CHKERR KSPSetDM(ksp, sub_dm);
CHKERR KSPSetDM(ksp, sub_dm);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSetUp(ksp);
auto L = createDMVector(sub_dm);
auto F = vectorDuplicate(L);
CHKERR KSPSolve(ksp, F, L);
CHKERR VecGhostUpdateBegin(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(L, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(sub_dm, L, INSERT_VALUES, SCATTER_REVERSE);
auto [error, th_error] = evaluateError();
MOFEM_LOG("LevelSet", Sev::inform) << "Error indicator " << error;
auto post_proc = [&](auto dm, auto out_name, auto th_error) {
MoFEMFunctionBegin;
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->setTagsToTransfer({th_error});
post_proc_fe->exeTestHook = test_mesh_bit;
if constexpr (DIM1 == 1 && DIM2 == 1) {
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
auto l_vec = boost::make_shared<VectorDouble>();
auto l_grad_mat = boost::make_shared<MatrixDouble>();
CHKERR AddHOOps<FE_DIM, FE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {L2});
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", l_vec));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("L", l_grad_mat));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"L", l_vec}},
{{"GradL", l_grad_mat}},
{}, {})
);
}
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
MoFEMFunctionReturn(0);
};
if constexpr (debug)
CHKERR post_proc(sub_dm, "dg_projection.h5m", th_error);
MoFEMFunctionReturn(0);
}
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
aggregate_projection_bit
constexpr int aggregate_projection_bit
all bits for projection problem
Definition: level_set.cpp:70
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
senseMap
std::array< int, 2 > senseMap
Definition: plate.cpp:102
Poisson2DiscontGalerkinOperators::indicesRowSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
Definition: PoissonDiscontinousGalerkin.hpp:23
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
A
constexpr AssemblyType A
Definition: level_set.cpp:32
LevelSet::readMesh
MoFEMErrorCode readMesh()
read mesh
Definition: level_set.cpp:502
MoFEM::CoreInterface::get_finite_element_meshset
virtual EntityHandle get_finite_element_meshset(const std::string name) const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMKSPSetComputeRHS
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
Definition: DMMoFEM.cpp:637
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
help
static char help[]
Definition: level_set.cpp:14
LevelSet::maxPtr
boost::shared_ptr< double > maxPtr
Definition: level_set.cpp:370
LevelSet::getSideFE
boost::shared_ptr< FaceSideEle > getSideFE(boost::shared_ptr< SideData > side_data_ptr)
create side element to assemble data from sides
Definition: level_set.cpp:997
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:887
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
NOISY
@ NOISY
Definition: definitions.h:224
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
LevelSet::OpMassLL
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< 1, DIM1 *DIM2 > OpMassLL
Definition: level_set.cpp:354
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
sdf_hertz_2d_axisymm_plane.zc
int zc
Definition: sdf_hertz_2d_axisymm_plane.py:8
LevelSet::OpRhsDomain::lPtr
boost::shared_ptr< MatrixDouble > lPtr
Definition: level_set.cpp:421
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition: DMMoFEM.hpp:1113
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
LevelSet::OpScalarFieldL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpBaseTimesVector< 1, DIM1 *DIM2, 1 > OpScalarFieldL
Definition: level_set.cpp:362
LevelSet::solveAdvection
MoFEMErrorCode solveAdvection()
solve advection problem
Definition: level_set.cpp:1934
LevelSet::OpRhsDomain::iNtegrate
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Definition: level_set.cpp:610
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
LevelSet::OpRhsSkeleton
Definition: level_set.cpp:435
FaceSideEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::FaceSideEle FaceSideEle
Definition: level_set.cpp:41
nb_levels
constexpr int nb_levels
Definition: level_set.cpp:58
LevelSet::OpLhsDomain::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
Definition: level_set.cpp:655
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
aggregate_bit
constexpr int aggregate_bit
all bits for advection problem
Definition: level_set.cpp:66
FaceSideEleOp
FaceSideEle::UserDataOperator FaceSideEleOp
Definition: level_set.cpp:45
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
LevelSet::initialiseFieldLevelSet
MoFEMErrorCode initialiseFieldLevelSet(boost::function< double(double, double, double)> level_fun=get_level_set)
initialise field set
Definition: level_set.cpp:1693
Poisson2DiscontGalerkinOperators::rowBaseSideMap
std::array< MatrixDouble, 2 > rowBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:26
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
LevelSet::OpSourceL
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< 1, DIM1 *DIM2 > OpSourceL
Definition: level_set.cpp:356
EntData
EntitiesFieldData::EntData EntData
Definition: level_set.cpp:36
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
potential_velocity_field_dim
constexpr size_t potential_velocity_field_dim
Definition: level_set.cpp:50
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEM::DMMoFEMKSPSetComputeOperators
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.
Definition: DMMoFEM.cpp:678
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
RIGHT_SIDE
@ RIGHT_SIDE
Definition: plate.cpp:92
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
I
FTensor::Index< 'I', DIM1 > I
Definition: level_set.cpp:29
FE_DIM
constexpr int FE_DIM
[Define dimension]
Definition: level_set.cpp:19
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
LevelSet::OpRhsDomain
Definition: level_set.cpp:412
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
LevelSet::testOp
MoFEMErrorCode testOp()
test consistency between tangent matrix and the right hand side vectors
Definition: level_set.cpp:1602
LevelSet::pushOpSkeleton
MoFEMErrorCode pushOpSkeleton()
push operator to integrate on skeleton
Definition: level_set.cpp:1173
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:251
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
SPACE_DIM
constexpr int SPACE_DIM
Definition: level_set.cpp:20
LevelSet::OpSourceV
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< G >::OpSource< potential_velocity_field_dim, potential_velocity_field_dim > OpSourceV
Definition: level_set.cpp:360
LevelSet::OpLhsDomain
Definition: level_set.cpp:426
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
current_bit
constexpr int current_bit
dofs bit used to do calculations
Definition: level_set.cpp:63
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
LevelSet::OpRhsSkeleton::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: level_set.cpp:703
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
LevelSet::get_level_set
static double get_level_set(const double x, const double y, const double z)
inital level set, i.e. advected filed
Definition: level_set.cpp:378
main
int main(int argc, char *argv[])
Definition: level_set.cpp:464
BoundaryEleOp
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
k
FTensor::Index< 'k', SPACE_DIM > k
Definition: level_set.cpp:581
MoFEM::DMSubDMSetUp_MoFEM
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
Definition: DMMoFEM.cpp:1350
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
DIM1
constexpr int DIM1
Definition: level_set.cpp:21
DomainEleOp
DomainEle::UserDataOperator DomainEleOp
Definition: level_set.cpp:43
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
LevelSet::setUpProblem
MoFEMErrorCode setUpProblem()
create fields, and set approximation order
Definition: level_set.cpp:560
LevelSet::OpLhsSkeleton::OpLhsSkeleton
OpLhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
Definition: level_set.cpp:604
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
LevelSet::OpLhsDomain::OpLhsDomain
OpLhsDomain(const std::string field_name, boost::shared_ptr< MatrixDouble > vel_ptr)
Definition: level_set.cpp:590
LevelSet::OpRhsDomain::velPtr
boost::shared_ptr< MatrixDouble > velPtr
Definition: level_set.cpp:423
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::TsCtx
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
start_bit
constexpr int start_bit
Definition: level_set.cpp:60
Poisson2DiscontGalerkinOperators::get_ntensor
auto get_ntensor(T &base_mat)
Definition: PoissonDiscontinousGalerkin.hpp:90
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
level_set_raw_ptr
LevelSet * level_set_raw_ptr
Definition: level_set.cpp:1931
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::PipelineManager::getOpDomainLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
Definition: PipelineManager.hpp:749
LevelSet::OpMassVV
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< G >::OpMass< potential_velocity_field_dim, potential_velocity_field_dim > OpMassVV
Definition: level_set.cpp:358
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
BiLinearForm
LevelSet::OpRhsSkeleton::OpRhsSkeleton
OpRhsSkeleton(boost::shared_ptr< SideData > side_data_ptr, boost::shared_ptr< FaceSideEle > side_fe_ptr)
Definition: level_set.cpp:598
MoFEM::OpCalculateTensor2FieldValuesDot
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:902
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index
Definition: Index.hpp:23
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::PipelineManager::getOpSkeletonLhsPipeline
boost::ptr_deque< UserDataOperator > & getOpSkeletonLhsPipeline()
Get the Op Skeleton Lhs Pipeline object.
Definition: PipelineManager.hpp:845
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: level_set.cpp:579
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:249
LevelSet::evaluateError
std::tuple< double, Tag > evaluateError()
evaluate error
Definition: level_set.cpp:1203
sdf_hertz_2d_axisymm_plane.yc
int yc
Definition: sdf_hertz_2d_axisymm_plane.py:7
LevelSet::runProblem
MoFEMErrorCode runProblem()
Definition: level_set.cpp:386
LevelSet::mField
MoFEM::Interface & mField
integrate skeleton operators on khs
Definition: level_set.cpp:346
Range
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
save_range
auto save_range
Definition: thermo_elastic.cpp:144
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FTensor::Tensor0
Definition: Tensor0.hpp:16
skeleton_bit
constexpr int skeleton_bit
skeleton elements bit
Definition: level_set.cpp:65
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
LevelSet::initialiseFieldVelocity
MoFEMErrorCode initialiseFieldVelocity(boost::function< double(double, double, double)> vel_fun=get_velocity_potential< FE_DIM >)
initialise potential velocity field
Definition: level_set.cpp:1814
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Poisson2DiscontGalerkinOperators::colBaseSideMap
std::array< MatrixDouble, 2 > colBaseSideMap
Definition: PoissonDiscontinousGalerkin.hpp:27
LevelSet::OpLhsSkeleton::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: level_set.cpp:809
LevelSet::refineMesh
MoFEMErrorCode refineMesh(WrapperClass &&wp)
Definition: level_set.cpp:2147
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
LevelSet::dgProjection
MoFEMErrorCode dgProjection(const int prj_bit=projection_bit)
dg level set projection
Definition: level_set.cpp:2401
debug
constexpr bool debug
Definition: level_set.cpp:53
G
constexpr IntegrationType G
Definition: level_set.cpp:33
MoFEM::DMMoFEMTSSetMonitor
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.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:1044
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Simple::getBitRefLevel
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:354
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
LEFT_SIDE
@ LEFT_SIDE
Definition: plate.cpp:92
LevelSet::OpRhsDomain::OpRhsDomain
OpRhsDomain(const std::string field_name, boost::shared_ptr< MatrixDouble > l_ptr, boost::shared_ptr< MatrixDouble > l_dot_ptr, boost::shared_ptr< MatrixDouble > vel_ptr)
Definition: level_set.cpp:583
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
DIM2
constexpr int DIM2
Definition: level_set.cpp:22
Poisson2DiscontGalerkinOperators::indicesColSideMap
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
Definition: PoissonDiscontinousGalerkin.hpp:25
LevelSet::LEFT_SIDE
@ LEFT_SIDE
Definition: level_set.cpp:367
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
LevelSet::getZeroLevelVelOp
ForcesAndSourcesCore::UserDataOperator * getZeroLevelVelOp(boost::shared_ptr< MatrixDouble > vel_ptr)
Get operator calculating velocity on coarse mesh.
Definition: level_set.cpp:922
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
LevelSet::pushOpDomain
MoFEMErrorCode pushOpDomain()
push operators to integrate operators on domain
Definition: level_set.cpp:954
sdf_hertz_2d_axisymm_plane.xc
int xc
Definition: sdf_hertz_2d_axisymm_plane.py:6
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
projection_bit
constexpr int projection_bit
Definition: level_set.cpp:68
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::make_array
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1990
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
j
FTensor::Index< 'j', SPACE_DIM > j
Definition: level_set.cpp:580
LevelSet
Definition: level_set.cpp:73
LevelSet::OpRhsDomain::lDotPtr
boost::shared_ptr< MatrixDouble > lDotPtr
Definition: level_set.cpp:422
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
ElementSide
ElementSide
Definition: plate.cpp:92
potential_velocity_space
constexpr FieldSpace potential_velocity_space
Definition: level_set.cpp:49
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
LevelSet::OpLhsSkeleton
Definition: level_set.cpp:450
LevelSet::RIGHT_SIDE
@ RIGHT_SIDE
Definition: level_set.cpp:367
MoFEM::OpCalculateHcurlVectorCurl
Calculate curl of vector field.
Definition: UserDataOperators.hpp:2488
DomianParentEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
Definition: level_set.cpp:39
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
LevelSet::testSideFE
MoFEMErrorCode testSideFE()
test integration side elements
Definition: level_set.cpp:1387
Generated by Doxygen 1.8.17 and hosted at University of Glasgow