v0.13.0
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
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Publicrule
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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 {
VectorDouble &vAls;
MatrixDouble &diffVals;
const bool checkGradients;
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
: DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
checkGradients(check_grads) {}
MoFEMErrorCode doWork(int side, EntityType type,
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;
}
if (checkGradients) {
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;
}
}
}
}
};
VectorDouble &vAls;
MatrixDouble &diffVals;
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);
}
MoFEMErrorCode doWork(int side, EntityType type,
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;
}
if (checkGradients) {
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(
MOFEM_LOG_TAG("AT", "atom_test");
// Create MoFEM instance
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
Simple *simple_interface = m_field.getInterface<Simple>();
PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
CHKERR simple_interface->getOptions();
CHKERR simple_interface->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_interface->addDomainField("FIELD1", space, base, 1);
CHKERR simple_interface->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_interface->setFieldOrder("FIELD1", approx_order + 1,
&rise_order);
CHKERR simple_interface->setFieldOrder("FIELD1", approx_order + 2,
&rise_twice);
}
CHKERR simple_interface->setUp();
auto dm = simple_interface->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 = [&]() {
if (SPACE_DIM == 2)
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOWeightsOnFace());
using OpSource = FormsIntegrators<DomainEleOp>::Assembly<
if (SPACE_DIM == 2) {
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOWeightsOnFace());
}
if (SPACE_DIM == 3) {
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJacVolume(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, nullptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOWeights(det_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacVolume(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOWeights(det_ptr));
}
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<
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_interface->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 = [&] {
boost::shared_ptr<VectorDouble> ptr_values(new VectorDouble());
boost::shared_ptr<MatrixDouble> ptr_diff_vals(new MatrixDouble());
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
if (SPACE_DIM == 2) {
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacSpaceForFaceImpl<2>(space, inv_jac_ptr));
}
if (SPACE_DIM == 3) {
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacVolume(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOInvJacToScalarBases(space, inv_jac_ptr));
}
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(
new OpCalculateScalarFieldGradient<SPACE_DIM>("FIELD1",
ptr_diff_vals));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
vals, diff_vals, ptr_values, ptr_diff_vals, true));
CHKERR pipeline_mng->loopFiniteElements();
};
CHKERR assemble_matrices_and_vectors();
CHKERR solve_problem();
CHKERR check_solution();
}
}
static Index< 'o', 3 > o
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
static int approx_order
static char help[]
constexpr int SPACE_DIM
static FTensor::Tensor1< double, 3 > fUn(const double x, const double y)
static FTensor::Tensor2< double, 3, 2 > diffFun(const double x, const double y)
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator