v0.13.1
scalar_check_approximation.cpp

Approximate polynomial in 2D and check derivatives

/**
* \file scalar_check_approximation.cpp
* \example scalar_check_approximation.cpp
*
* Approximate polynomial in 2D and check derivatives
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
static int approx_order = 4;
template <int DIM> struct ApproxFunctionsImpl {};
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
template <> struct ApproxFunctionsImpl<2> {
static double fUn(const double x, const double y, double z) {
double r = 1;
for (int o = 1; o <= approx_order; ++o) {
for (int i = 0; i <= o; ++i) {
int j = o - i;
if (j >= 0)
r += pow(x, i) * pow(y, j);
}
}
return r;
}
static FTensor::Tensor1<double, 2> diffFun(const double x, const double y,
double z) {
for (int o = 1; o <= approx_order; ++o) {
for (int i = 0; i <= o; ++i) {
int j = o - i;
if (j >= 0) {
r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) : 0;
r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) : 0;
}
}
}
return r;
}
};
template <> struct ApproxFunctionsImpl<3> {
static double fUn(const double x, const double y, double z) {
double r = 1;
for (int o = 1; o <= approx_order; ++o) {
for (int i = 0; i <= o; ++i) {
for (int j = 0; j <= o - i; j++) {
int k = o - i - j;
if (k >= 0) {
r += pow(x, i) * pow(y, j) * pow(z, k);
}
}
}
}
return r;
}
static FTensor::Tensor1<double, 3> diffFun(const double x, const double y,
double z) {
for (int o = 1; o <= approx_order; ++o) {
for (int i = 0; i <= o; ++i) {
for (int j = 0; j <= o - i; j++) {
int k = o - i - j;
if (k >= 0) {
r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) * pow(z, k) : 0;
r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) * pow(z, k) : 0;
r(2) += k > 0 ? k * pow(x, i) * pow(y, j) * pow(z, k - 1) : 0;
}
}
}
}
return r;
}
};
struct OpValsDiffVals : public DomainEleOp {
const bool checkGradients;
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
: DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
checkGradients(check_grads) {}
const int nb_gauss_pts = getGaussPts().size2();
if (type == MBVERTEX) {
vAls.resize(nb_gauss_pts, false);
diffVals.resize(SPACE_DIM, nb_gauss_pts, false);
vAls.clear();
diffVals.clear();
}
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
MOFEM_LOG("AT", Sev::noisy)
<< "Type " << moab::CN::EntityTypeName(type) << " side " << side;
MOFEM_LOG("AT", Sev::noisy) << data.getN();
MOFEM_LOG("AT", Sev::noisy) << data.getDiffN();
auto t_vals = getFTensor0FromVec(vAls);
auto t_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int bb = 0; bb != nb_dofs; bb++) {
t_vals += t_base_fun * t_data;
++t_base_fun;
++t_data;
}
const double v = t_vals;
if (!std::isnormal(v))
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
++t_vals;
}
auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
auto t_diff_base_fun = data.getFTensor1DiffN<SPACE_DIM>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int bb = 0; bb != nb_dofs; bb++) {
t_diff_vals(i) += t_diff_base_fun(i) * t_data;
++t_diff_base_fun;
++t_data;
}
for (int d = 0; d != SPACE_DIM; ++d)
if (!std::isnormal(t_diff_vals(d)))
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
++t_diff_vals;
}
}
}
}
};
boost::shared_ptr<VectorDouble> ptrVals;
boost::shared_ptr<MatrixDouble> ptrDiffVals;
const bool checkGradients;
boost::shared_ptr<VectorDouble> &ptr_vals,
boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
bool check_grads)
: DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
checkGradients(check_grads) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
const double eps = 1e-6;
const int nb_gauss_pts = data.getN().size1();
auto t_vals = getFTensor0FromVec(vAls);
auto t_ptr_vals = getFTensor0FromVec(*ptrVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double err_val;
// Check user data operators
err_val = std::abs(t_vals - t_ptr_vals);
MOFEM_LOG("AT", Sev::noisy) << "Val op error " << err_val;
if (err_val > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong value from operator %4.3e", err_val);
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
const double z = getCoordsAtGaussPts()(gg, 2);
// Check approximation
const double delta_val = t_vals - ApproxFunctions::fUn(x, y, z);
err_val = std::fabs(delta_val * delta_val);
MOFEM_LOG("AT", Sev::verbose) << err_val << " : " << t_vals;
if (err_val > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong value %4.3e",
err_val);
++t_vals;
++t_ptr_vals;
}
auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double err_diff_val;
t_delta_diff_val(i) = t_diff_vals(i) - t_ptr_diff_vals(i);
err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
MOFEM_LOG("AT", Sev::noisy) << "Diff val op error " << err_diff_val;
if (err_diff_val > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong derivatives from operator %4.3e", err_diff_val);
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
const double z = getCoordsAtGaussPts()(gg, 2);
// Check approximation
auto t_diff_anal = ApproxFunctions::diffFun(x, y, z);
t_delta_diff_val(i) = t_diff_vals(i) - t_diff_anal(i);
err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
if (SPACE_DIM == 3)
MOFEM_LOG("AT", Sev::noisy)
<< "Diff val " << err_diff_val << " : "
<< sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
<< t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
<< t_diff_vals(1) << " (" << t_diff_anal(1) << ") "
<< t_diff_vals(2) << " (" << t_diff_anal(2) << ")";
else
MOFEM_LOG("AT", Sev::noisy)
<< "Diff val " << err_diff_val << " : "
<< sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
<< t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
<< t_diff_vals(1) << " (" << t_diff_anal(1) << ")";
MOFEM_LOG("AT", Sev::verbose)
<< getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
MOFEM_LOG("AT", Sev::verbose)
<< getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
MOFEM_LOG("AT", Sev::verbose)
<< getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
MOFEM_LOG("AT", Sev::verbose) << "Diff val error " << err_diff_val;
if (err_diff_val > eps)
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong derivative of value %4.3e %4.3e", err_diff_val,
t_diff_anal.l2());
++t_diff_vals;
++t_ptr_diff_vals;
}
}
}
};
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "AT"));
LogManager::setLog("AT");
MOFEM_LOG_TAG("AT", "atom_test");
// Create MoFEM instance
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
CHKERR simple->getOptions();
simple->getAddBoundaryFE() = true;
CHKERR simple->loadFile("", "");
// Declare elements
enum bases {
AINSWORTH,
AINSWORTH_LOBATTO,
DEMKOWICZ,
BERNSTEIN,
LASBASETOP
};
const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
if (choice_base_value == AINSWORTH)
if (choice_base_value == AINSWORTH_LOBATTO)
else if (choice_base_value == DEMKOWICZ)
else if (choice_base_value == BERNSTEIN)
enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
const char *list_spaces[] = {"h1", "l2"};
PetscInt choice_space_value = H1SPACE;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
LASBASETSPACE, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
FieldSpace space = H1;
if (choice_space_value == H1SPACE)
space = H1;
else if (choice_space_value == L2SPACE)
space = L2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
PETSC_NULL);
CHKERR simple->addDomainField("FIELD1", space, base, 1);
CHKERR simple->setFieldOrder("FIELD1", approx_order);
Range edges, faces;
CHKERR moab.get_entities_by_dimension(0, 1, edges);
CHKERR moab.get_entities_by_dimension(0, 2, faces);
if (choice_base_value != BERNSTEIN) {
Range rise_order;
int i = 0;
for (auto e : edges) {
if (!(i % 2)) {
rise_order.insert(e);
}
++i;
}
for (auto f : faces) {
if (!(i % 3)) {
rise_order.insert(f);
}
++i;
}
Range rise_twice;
for (auto e : rise_order) {
if (!(i % 2)) {
rise_twice.insert(e);
}
++i;
}
CHKERR simple->setFieldOrder("FIELD1", approx_order + 1, &rise_order);
CHKERR simple->setFieldOrder("FIELD1", approx_order + 2, &rise_twice);
}
CHKERR simple->defineFiniteElements();
auto volume_adj = [](moab::Interface &moab, const Field &field,
const EntFiniteElement &fe,
std::vector<EntityHandle> &adjacency) {
EntityHandle fe_ent = fe.getEnt();
switch (field.getSpace()) {
case H1:
CHKERR moab.get_connectivity(&fe_ent, 1, adjacency, true);
case HCURL:
CHKERR moab.get_adjacencies(&fe_ent, 1, 1, false, adjacency,
moab::Interface::UNION);
case HDIV:
case L2:
CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, adjacency,
moab::Interface::UNION);
adjacency.push_back(fe_ent);
// build side table
for (auto ent : adjacency)
fe.getSideNumberPtr(ent);
break;
case NOFIELD: {
CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
false);
for (auto ent : adjacency) {
const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
.insert(
boost::shared_ptr<SideNumber>(new SideNumber(ent, -1, 0, 0)));
}
} break;
default:
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"this field is not implemented for TRI finite element");
}
};
simple->getDomainFEName(), MBTET, volume_adj);
simple->getDomainFEName(), MBHEX, volume_adj);
CHKERR simple->defineProblem(PETSC_TRUE);
CHKERR simple->buildFields();
CHKERR simple->buildFiniteElements();
CHKERR simple->buildProblem();
auto dm = simple->getDM();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
MatrixDouble diff_vals;
auto assemble_matrices_and_vectors = [&]() {
pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
"FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSource("FIELD1", ApproxFunctions::fUn));
auto integration_rule = [](int, int, int p_data) {
return 2 * p_data + 1;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
};
auto solve_problem = [&] {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
};
auto check_solution = [&] {
auto ptr_values = boost::make_shared<VectorDouble>();
auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline(), {space});
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpValsDiffVals(vals, diff_vals, true));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldValues("FIELD1", ptr_values));
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_diff_vals));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
vals, diff_vals, ptr_values, ptr_diff_vals, true));
CHKERR pipeline_mng->loopFiniteElements();
};
auto post_proc = [&] {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
post_proc_fe->getOpPtrVector(), {space});
auto ptr_values = boost::make_shared<VectorDouble>();
auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("FIELD1", ptr_values));
post_proc_fe->getOpPtrVector().push_back(
ptr_diff_vals));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"FIELD1", ptr_values}},
{{"FIELD1_GRAD", ptr_diff_vals}},
{}, {})
);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile("out_post_proc.h5m");
auto bdy_post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<BoundaryEle>>(m_field);
auto op_loop_side = new OpLoopSide<EleOnSide>(
m_field, simple->getDomainFEName(), SPACE_DIM);
// push operators to side element
op_loop_side->getOpPtrVector(), {space});
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("FIELD1", ptr_values));
op_loop_side->getOpPtrVector().push_back(
ptr_diff_vals));
// push op to boundary element
bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
bdy_post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
bdy_post_proc_fe->getPostProcMesh(),
bdy_post_proc_fe->getMapGaussPts(),
{{"FIELD1", ptr_values}},
{{"FIELD1_GRAD", ptr_diff_vals}},
{}, {})
);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
bdy_post_proc_fe);
bdy_post_proc_fe->writeFile("out_post_proc_bdy.h5m");
};
CHKERR assemble_matrices_and_vectors();
CHKERR solve_problem();
CHKERR check_solution();
CHKERR post_proc();
}
}
static Index< 'o', 3 > o
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
const double r
rate factor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
static int approx_order
static char help[]
ElementsAndOps< SPACE_DIM >::EleOnSide EleOnSide
constexpr int SPACE_DIM
static FTensor::Tensor2< double, BASE_DIM, SPACE_DIM > diffFun(const double x, const double y)
static FTensor::Tensor1< double, BASE_DIM > fUn(const double x, const double y, double z)
Add operators pushing bases from local to physical configuration.
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Finite element data for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Base face element used to integrate on skeleton.
Provide data structure for (tensor) field approximation.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Base volume element used to integrate on skeleton.
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > ptrDiffVals
OpCheckValsDiffVals(boost::shared_ptr< MatrixDouble > ptr_vals, boost::shared_ptr< VectorDouble > ptr_div, boost::shared_ptr< MatrixDouble > ptr_grad, boost::shared_ptr< MatrixDouble > ptr_hess)
boost::shared_ptr< MatrixDouble > ptrVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
FTensor::Index< 'i', SPACE_DIM > i