Testing Discontinuous Galerkin (DG) projection operators.
static char help[] =
"DG Projection Test - validates discontinuous Galerkin "
"projection accuracy\n\n";
auto fun = [](
const double x,
const double y,
const double z) {
return x + y + x * x + y * y;
};
private:
};
};
auto save_range = [](moab::Interface &moab,
const std::string name,
const Range r, std::vector<Tag> tags = {}) {
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1,
tags.data(), tags.size());
} else {
MOFEM_LOG(
"SELF", Sev::warning) <<
"Empty range for " << name;
}
};
OpError(boost::shared_ptr<MatrixDouble> data_ptr)
const int nb_integration_pts = getGaussPts().size2();
auto t_val = getFTensor1FromMat<1>(*(
dataPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double projected_value = t_val(0);
double analytical_value =
fun(t_coords(0), t_coords(1), t_coords(2));
double error = projected_value - analytical_value;
constexpr double eps = 1e-8;
if (std::abs(error) >
eps) {
<< "Projection error too large: " << error << " at point ("
<< t_coords(0) << ", " << t_coords(1) << ")"
<< " projected=" << projected_value
<< " analytical=" << analytical_value;
"DG projection failed accuracy test");
}
++t_val;
++t_coords;
}
MOFEM_LOG(
"SELF", Sev::noisy) <<
"DG projection accuracy validation passed";
}
private:
boost::shared_ptr<MatrixDouble>
dataPtr;
};
}
char mesh_File_Name[255];
mesh_File_Name, 255, PETSC_NULLPTR);
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solving DG projection system";
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
auto rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
auto get_child_op = [&](auto &pip) {
mField,
simple->getDomainFEName(), child_bit, child_bit, Sev::noisy);
auto fe_child_ptr = op_this_child->getThisFEPtr();
fe_child_ptr->getRuleHook = [] (
int,
int,
int p) {
return -1; };
child_bit, child_bit, MBEDGE, child_edges);
CHKERR setDGSetIntegrationPoints<SPACE_DIM>(
fe_child_ptr->setRuleHook, [](int, int, int p) { return 2 * p; },
boost::make_shared<Range>(child_edges));
pip.push_back(op_this_child);
return fe_child_ptr;
};
auto get_filed_eval_op = [&](auto fe_child_ptr) {
parent_bit);
auto data_ptr = boost::make_shared<MatrixDouble>();
auto eval_data_ptr = boost::make_shared<MatrixDouble>();
if (auto fe_eval_ptr = field_eval_data->feMethodPtr.lock()) {
fe_eval_ptr->getRuleHook = [] (
int,
int,
int p) {
return -1; };
fe_eval_ptr->getOpPtrVector().push_back(
eval_data_ptr));
op_test->doWorkRhsHook =
<< "Field evaluator method pointer is valid";
<< op_ptr->getGaussPts();
<< "Loop size " << op_ptr->getPtrFE()->getLoopSize();
<< "Coords at gauss pts: " << op_ptr->getCoordsAtGaussPts();
};
fe_eval_ptr->getOpPtrVector().push_back(op_test);
} else {
"Field evaluator method pointer is expired");
}
auto op_ptr = field_eval_ptr->getDataOperator<
SPACE_DIM>(
fe_child_ptr->getOpPtrVector().push_back(op_ptr);
return data_ptr;
};
auto dg_projection_base = [&](auto fe_child_ptr, auto data_ptr) {
constexpr int projection_order =
order;
auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto mass_ptr = boost::make_shared<MatrixDouble>();
auto coeffs_ptr = boost::make_shared<MatrixDouble>();
fe_child_ptr->getOpPtrVector().push_back(
fe_child_ptr->getOpPtrVector().push_back(
new OpError(data_ptr));
op_set_coeffs->doWorkRhsHook =
if (!field_ents.size())
if (auto e_ptr = field_ents[0]) {
auto field_ent_data = e_ptr->getEntFieldData();
std::copy(coeffs_ptr->data().begin(), coeffs_ptr->data().end(),
field_ent_data.begin());
}
};
fe_child_ptr->getOpPtrVector().push_back(op_set_coeffs);
};
auto fe_child_ptr = get_child_op(pipeline_mng->getOpDomainRhsPipeline());
CHKERR dg_projection_base(fe_child_ptr, get_filed_eval_op(fe_child_ptr));
auto test_data_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpError(test_data_ptr));
CHKERR pipeline_mng->loopFiniteElements();
}
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceEle>(
mField);
post_proc_fe->getOpPtrVector(), {H1});
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{FIELD_NAME, u_ptr}},
OpPPMap::DataMapMat{{
"GRAD_" + std::string(
FIELD_NAME), grad_u_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(file_name);
}
auto make_edge_flip = [&](
auto edge,
auto adj_faces,
Range &new_tris) {
int num_nodes;
CHKERR moab.get_connectivity(e, conn, num_nodes,
true);
std::copy(conn, conn + num_nodes, conn_cpy);
};
std::array<EntityHandle, 6> adj_conn;
CHKERR get_conn(adj_faces[0], &adj_conn[0]);
CHKERR get_conn(adj_faces[1], &adj_conn[3]);
std::array<EntityHandle, 2> edge_conn;
CHKERR get_conn(edge, edge_conn.data());
std::array<EntityHandle, 2> new_edge_conn;
for (
int i = 0;
i != 6; ++
i) {
if (adj_conn[
i] != edge_conn[0] && adj_conn[
i] != edge_conn[1]) {
new_edge_conn[
j] = adj_conn[
i];
}
}
auto &new_conn = adj_conn;
for (
int t = 0;
t != 2; ++
t) {
for (
int i = 0;
i != 3; ++
i) {
if (
(adj_conn[3 *
t +
i % 3] == edge_conn[0] &&
adj_conn[3 *
t + (
i + 1) % 3] == edge_conn[1])
||
(adj_conn[3 *
t +
i % 3] == edge_conn[1] &&
adj_conn[3 *
t + (
i + 1) % 3] == edge_conn[0])
) {
new_conn[3 *
t + (
i + 1) % 3] = new_edge_conn[
t];
break;
}
}
}
for (
int t = 0;
t != 2; ++
t) {
false, rtri);
if (!rtri.size()) {
new_tris.insert(tri);
} else {
#ifndef NDEBUG
if (rtri.size() != 1) {
<< "Multiple tries created during edge flip for edge " << edge
<< " adjacent faces " << std::endl
<< rtri;
"Multiple tries created during edge flip");
}
#endif
new_tris.merge(rtri);
}
}
moab::Interface::UNION);
#ifndef NDEBUG
#endif
};
Skinner skin(&moab);
CHKERR skin.find_skin(0, tris,
false, skin_edges);
Range new_tris, flipped_tris;
int flip_count = 0;
for (auto edge : edges) {
adjacent_tris = subtract(adjacent_tris, adj_skin_tris);
adjacent_tris = subtract(adjacent_tris, flipped_tris);
if (adjacent_tris.size() == 2) {
CHKERR make_edge_flip(edge, adjacent_tris, new_tris);
flipped_tris.merge(adjacent_tris);
++flip_count;
}
}
Range not_flipped_tris = subtract(all_tris, flipped_tris);
<< "Flipped " << flip_count << " edges with two adjacent faces.";
child_bit);
child_bit);
child_bit,
BitRefLevel().set(),
"edge_flips_before_refinement.vtk",
"VTK",
"");
}
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
Order displacement.
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFaceEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
constexpr auto field_name
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxGradVals
boost::shared_ptr< MatrixDouble > approxHessianVals
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode projectResults(BitRefLevel parent_bit, BitRefLevel child_bit)
[Solve]
MoFEMErrorCode reSetupProblem(BitRefLevel child_bit)
[Edge flips]
MoFEMErrorCode edgeFlips(BitRefLevel parent_bit, BitRefLevel child_bit)
[Output results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
Reference to MoFEM interface.
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
Field evaluator interface.
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Specialization for double precision vector field values calculation.
Execute "this" element in the operator.
Post post-proc data at points from hash maps.
Template struct for dimension-specific finite element types.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain left-hand side finite element.
Simple interface for fast problem set-up.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode reSetUp(bool only_dm=false)
Rebuild internal MoFEM data structures.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
DomainEle::UserDataOperator DomainEleOp