v0.13.1
hcurl_check_approx_in_2d.cpp

Approximate polynomial in 2D and check derivatives

/**
* \file hcurl_check_approx_in_2d
* \example hcurl_check_approx_in_2d.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 Public
* 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";
constexpr int BASE_DIM = 3;
constexpr int SPACE_DIM = 2;
/**
* @brief OPerator to integrate mass matrix for least square approximation
*
*/
/**
* @brief Operator to integrate the right hand side matrix for the problem
*
*/
GAUSS>::OpSource<BASE_DIM, BASE_DIM>;
constexpr double a0 = 0.0;
constexpr double a1 = 2.0;
constexpr double a2 = -15.0 * a0;
constexpr double a3 = -20.0 / 6 * a1;
constexpr double a4 = 15.0 * a0;
constexpr double a5 = a1;
constexpr double a6 = -a0;
static FTensor::Tensor1<double, BASE_DIM> fUn(const double x, const double y,
double z) {
// x
6 * a6 * std::pow(x, 5) * std::pow(y, 0) +
5 * a5 * std::pow(x, 4) * std::pow(y, 1) +
4 * a4 * std::pow(x, 3) * std::pow(y, 2) +
3 * a3 * std::pow(x, 2) * std::pow(y, 3) +
2 * a2 * std::pow(x, 1) * std::pow(y, 4) +
1 * a1 * std::pow(x, 0) * std::pow(y, 5),
// y
1 * a5 * std::pow(x, 5) * std::pow(y, 0) +
2 * a4 * std::pow(x, 4) * std::pow(y, 1) +
3 * a3 * std::pow(x, 3) * std::pow(y, 2) +
4 * a2 * std::pow(x, 2) * std::pow(y, 3) +
5 * a1 * std::pow(x, 1) * std::pow(y, 4) +
6 * a0 * std::pow(x, 0) * std::pow(y, 5),
// z
0.);
}
const double y) {
// x,x
30 * a6 * pow(x, 4) * pow(y, 0) + 20 * a5 * pow(x, 3) * pow(y, 1) +
12 * a4 * pow(x, 2) * pow(y, 2) + 6 * a3 * pow(x, 1) * pow(y, 3) +
2 * a2 * pow(x, 0) * pow(y, 4),
// x,y
5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
5 * a1 * pow(x, 0) * pow(y, 4),
// y,x
5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
5 * a1 * pow(x, 0) * pow(y, 4),
// y,y
2 * a4 * pow(x, 4) * pow(y, 0) + 6 * a3 * pow(x, 3) * pow(y, 1) +
12 * a2 * pow(x, 2) * pow(y, 2) + 20 * a1 * pow(x, 1) * pow(y, 3) +
30 * a0 * pow(x, 0) * pow(y, 4),
// z
0., 0.);
}
diff2Fun(const double x, const double y) {
// x,xx 0/000
30 * 4 * a6 * pow(x, 3) * pow(y, 0) +
20 * 3 * a5 * pow(x, 2) * pow(y, 1) +
12 * 2 * a4 * pow(x, 1) * pow(y, 2) +
6 * 1 * a3 * pow(x, 0) * pow(y, 3),
// x,xy 1/001
20 * 1 * a5 * pow(x, 3) * pow(y, 0) +
12 * 2 * a4 * pow(x, 2) * pow(y, 2) +
6 * 3 * a3 * pow(x, 1) * pow(y, 2) +
2 * 4 * a2 * pow(x, 0) * pow(y, 3),
// x,yx 2/010
5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
8 * 1 * a2 * pow(x, 0) * pow(y, 3),
// x,yy 3/011
8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
5 * 4 * a1 * pow(x, 0) * pow(y, 3),
// y,xx 4/100
5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
8 * 1 * a2 * pow(x, 0) * pow(y, 3),
// y,xy 5/101
8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
5 * 4 * a1 * pow(x, 0) * pow(y, 3),
// y,yx 6/110
2 * 4 * a4 * pow(x, 3) * pow(y, 0) +
6 * 3 * a3 * pow(x, 2) * pow(y, 1) +
12 * 2 * a2 * pow(x, 1) * pow(y, 2) +
20 * 1 * a1 * pow(x, 0) * pow(y, 3),
// y,yy 7/111
6 * 1 * a3 * pow(x, 3) * pow(y, 0) +
12 * 2 * a2 * pow(x, 2) * pow(y, 1) +
20 * 3 * a1 * pow(x, 1) * pow(y, 2) +
30 * 4 * a0 * pow(x, 0) * pow(y, 3),
// z,xx 8/200
0.,
// z,xy 9/201
0.,
// z,yx 10/210
0.,
// z,yy 11/211
0.);
}
};
struct OpCheckValsDiffVals : public FaceEleOp {
boost::shared_ptr<MatrixDouble> ptrVals;
boost::shared_ptr<VectorDouble> ptrDiv;
boost::shared_ptr<MatrixDouble> ptrGrad;
boost::shared_ptr<MatrixDouble> ptrHess;
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)
: FaceEleOp("FIELD1", OPROW), ptrVals(ptr_vals), ptrDiv(ptr_div),
ptrGrad(ptr_grad), ptrHess(ptr_hess) {}
const double eps = 1e-6;
if (type == MBEDGE && side == 0) {
const int nb_gauss_pts = data.getN().size1();
auto t_vals_from_op = getFTensor1FromMat<3>(*ptrVals);
auto t_div_from_op = getFTensor0FromVec(*ptrDiv);
auto t_grad_from_op = getFTensor2FromMat<3, 2>(*ptrGrad);
auto t_hess_from_op = getFTensor3FromMat<3, 2, 2>(*ptrHess);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
// Check approximation
delta_val(i) = t_vals_from_op(i) - ApproxFunctions::fUn(x, y, 0)(i);
double err_val = sqrt(delta_val(i) * delta_val(i));
if (err_val > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong value %4.3e", err_val);
delta_diff_val(i, j) =
t_grad_from_op(i, j) - ApproxFunctions::diffFun(x, y)(i, j);
double err_diff_val = sqrt(delta_diff_val(i, j) * delta_diff_val(i, j));
if (err_diff_val > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong derivative of value %4.3e", err_diff_val);
double div = t_grad_from_op(0, 0) + t_grad_from_op(1, 1);
double err_div = div - t_div_from_op;
if (err_div > eps)
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
err_div, div, t_div_from_op);
delta_diff2_val(i, j, k) =
t_hess_from_op(i, j, k) - ApproxFunctions::diff2Fun(x, y)(i, j, k);
double hess_diff_error =
sqrt(delta_diff2_val(i, j, k) * delta_diff2_val(i, j, k));
if (hess_diff_error > eps)
SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"Wrong hessian from operator %4.3e", hess_diff_error);
++t_vals_from_op;
++t_div_from_op;
++t_grad_from_op;
++t_hess_from_op;
}
}
}
};
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;
// 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("", "rectangle_tri.h5m");
// Declare elements
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
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)
else if (choice_base_value == DEMKOWICZ)
CHKERR simple_interface->addDomainField("FIELD1", HCURL, base, 1);
constexpr int order = 5;
CHKERR simple_interface->setFieldOrder("FIELD1", order);
CHKERR simple_interface->setUp();
auto dm = simple_interface->getDM();
MatrixDouble vals, diff_vals;
auto assemble_matrices_and_vectors = [&]() {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("FIELD1", "FIELD1",
[](double, double, double) { return 1; })
);
auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
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 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 = [&] {
pipeline_mng->getOpDomainLhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().clear();
auto ptr_values = boost::make_shared<MatrixDouble>();
auto ptr_divergence = boost::make_shared<VectorDouble>();
auto ptr_grad = boost::make_shared<MatrixDouble>();
auto ptr_hessian = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
// Change H-curl to H-div in 2D, and apply Piola transform
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacHcurlFace(inv_jac_ptr));
// Evaluate base function second derivative
auto base_mass = boost::make_shared<MatrixDouble>();
auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2, base, L2));
pipeline_mng->getOpDomainRhsPipeline().push_back(
inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
base_mass, data_l2, base, HCURL));
// Calculate field values at integration points
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHVecVectorField<BASE_DIM>("FIELD1", ptr_values));
// Gradient
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_grad));
// Hessian
pipeline_mng->getOpDomainRhsPipeline().push_back(
"FIELD1", ptr_divergence));
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_hessian));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
CHKERR pipeline_mng->loopFiniteElements();
};
CHKERR assemble_matrices_and_vectors();
CHKERR solve_problem();
CHKERR check_solution();
}
}
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#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
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#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
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:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
constexpr double a5
constexpr double a2
int main(int argc, char *argv[])
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
static char help[]
constexpr double a4
constexpr int SPACE_DIM
constexpr double a0
constexpr int BASE_DIM
constexpr double a3
constexpr double a1
constexpr double a6
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
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
const double D
diffusivity
static FTensor::Tensor2< double, BASE_DIM, SPACE_DIM > diffFun(const double x, const double y)
static FTensor::Tensor3< double, BASE_DIM, SPACE_DIM, SPACE_DIM > diff2Fun(const double x, const double y)
static FTensor::Tensor1< double, BASE_DIM > fUn(const double x, const double y, double z)
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)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
Get vector field for H-div approximation.
Calculate gradient of vector field.
Calculate gradient of vector field.
Calculate divergence of vector field.
Make Hdiv space from Hcurl space in 2d.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
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.
Definition: Simple.cpp:378
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:294
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:810
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:308
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:593
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:753
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > ptrHess
boost::shared_ptr< MatrixDouble > ptrGrad
boost::shared_ptr< VectorDouble > ptrDiv
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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)
FTensor::Index< 'j', 2 > j
boost::shared_ptr< MatrixDouble > ptrVals
FTensor::Index< 'k', 2 > k