#include <BasicFiniteElements.hpp>
};
-1;
auto u_exact = [](
const double x,
const double y,
const double) {
return x * x * y * y;
else
return cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
else
-2 * M_PI * cos(2 * M_PI * y) * sin(2 * M_PI * x),
-2 * M_PI * cos(2 * M_PI * x) * sin(2 * M_PI * y)
};
};
auto source = [](
const double x,
const double y,
const double) {
return -(2 * x * x + 2 * y * y);
else
return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
static char help[] =
"...\n\n";
public:
private:
};
: domainField("U"), mField(m_field), oRder(4) {}
}
PETSC_NULL);
PETSC_NULL);
auto save_shared = [&](auto meshset, std::string prefix) {
auto file_name =
prefix + "_" +
1);
};
}
}
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
[](const double, const double, const double) { return 1; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
pipeline_mng->getOpSkeletonLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
auto rule_lhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto rule_rhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
}
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto u_vals_ptr = boost::make_shared<VectorDouble>();
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
enum NORMS {
L2 = 0, SEMI_NORM,
H1, LAST_NORM };
std::array<int, LAST_NORM> error_indices;
for (
int i = 0;
i != LAST_NORM; ++
i)
if (const size_t nb_dofs = data.getIndices().size()) {
const int nb_integration_pts =
o->getGaussPts().size2();
auto t_w =
o->getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
auto t_coords =
o->getFTensor1CoordsAtGaussPts();
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w *
o->getMeasure();
const double diff =
t_val -
u_exact(t_coords(0), t_coords(1), t_coords(2));
auto t_exact_grad =
u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
const double diff_grad =
(t_grad(
i) - t_exact_grad(
i)) * (t_grad(
i) - t_exact_grad(
i));
error[SEMI_NORM] +=
alpha * diff_grad;
++t_w;
++t_val;
++t_grad;
++t_coords;
}
error[
H1] = error[
L2] + error[SEMI_NORM];
ADD_VALUES);
}
};
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
std::array<VectorDouble, 2> side_vals;
std::array<double, 2> area_map;
side_vals_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
const auto nb_in_loop =
o->getFEMethod()->nInTheLoop;
area_map[nb_in_loop] =
o->getMeasure();
side_vals[nb_in_loop] = *u_vals_ptr;
if (!nb_in_loop) {
area_map[1] = 0;
side_vals[1].clear();
}
};
side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
CHKERR o->loopSideFaces(
"dFE", *side_fe_ptr);
const auto in_the_loop = side_fe_ptr->nInTheLoop;
#ifndef NDEBUG
const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
<< "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
#endif
const double s =
o->getMeasure() / (area_map[0] + area_map[1]);
constexpr std::array<int, 2> sign_array{1, -1};
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
const int nb_integration_pts =
o->getGaussPts().size2();
if (!in_the_loop) {
side_vals[1].resize(nb_integration_pts, false);
auto t_coords =
o->getFTensor1CoordsAtGaussPts();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
t_val_m =
u_exact(t_coords(0), t_coords(1), t_coords(2));
++t_coords;
++t_val_m;
}
}
auto t_w =
o->getFTensor0IntegrationWeight();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha =
o->getMeasure() * t_w;
const auto diff = t_val_p - t_val_m;
error[SEMI_NORM] +=
alpha *
p * diff * diff;
++t_w;
++t_val_p;
++t_val_m;
}
error[
H1] = error[SEMI_NORM];
ADD_VALUES);
};
skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
boundary_error_op->doWorkRhsHook = do_work_rhs_error;
pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(l2_vec);
CHKERR VecAssemblyEnd(l2_vec);
const double *array;
CHKERR VecGetArrayRead(l2_vec, &array);
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm L2 %6.4e",
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm Energetic %6.4e",
std::sqrt(array[SEMI_NORM]));
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm H1 %6.4e",
constexpr double eps = 1e-12;
if (std::sqrt(array[
H1]) >
eps)
}
CHKERR VecRestoreArrayRead(l2_vec, &array);
}
}
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
post_proc_fe->generateReferenceElementMesh();
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(
"out_result.h5m");
}
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
CHKERR poisson_problem.runProgram();
}
return 0;
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
@ 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.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
constexpr auto VecSetValues
const double D
diffusivity
int main(int argc, char *argv[])
[Run program]
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEleOp DomainEleOp
ElementsAndOps< SPACE_DIM >::BoundaryEleOp BoundaryEleOp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainGradGrad
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
Base face element used to integrate on skeleton.
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
keeps basic data about problem
DofIdx getNbDofsRow() const
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
MoFEM::Interface & mField
MoFEMErrorCode boundaryCondition()
[Setup problem]
MoFEMErrorCode checkResults()
[Solve system]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode outputResults()
[Output results]
MoFEMErrorCode runProgram()
[Output results]
Opator tp evaluate Dirichlet boundary conditions using DG.
Operator the left hand side matrix.