v0.13.1
Loading...
Searching...
No Matches
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
*
*/
#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(
new OpDomainSource("FIELD1", ApproxFunctions::fUn));
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 char help[]
int main()
Definition: adol-c_atom.cpp:46
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int SPACE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#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
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#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
#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
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
constexpr double a5
constexpr double a2
constexpr double a4
constexpr double a0
constexpr int BASE_DIM
constexpr double a3
constexpr double a1
constexpr double a6
double D
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
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: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
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:26
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:373
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:289
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:780
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:303
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:588
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:723
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.
FTensor::Index< 'j', 2 > j
boost::shared_ptr< MatrixDouble > ptrVals
FTensor::Index< 'k', 2 > k