static char help[] =
"...\n\n";
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;
double z) {
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),
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),
0.);
}
const double y) {
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),
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),
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),
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),
0., 0.);
}
diff2Fun(
const double x,
const double y) {
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),
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),
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),
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),
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),
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),
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),
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),
0.,
0.,
0.,
0.);
}
};
boost::shared_ptr<MatrixDouble>
ptrVals;
boost::shared_ptr<VectorDouble>
ptrDiv;
boost::shared_ptr<MatrixDouble>
ptrGrad;
boost::shared_ptr<MatrixDouble>
ptrHess;
boost::shared_ptr<VectorDouble> ptr_div,
boost::shared_ptr<MatrixDouble> ptr_grad,
boost::shared_ptr<MatrixDouble> ptr_hess)
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++) {
double err_val = sqrt(delta_val(
i) * delta_val(
i));
"Wrong value %4.3e", err_val);
double err_diff_val = sqrt(delta_diff_val(
i,
j) * delta_diff_val(
i,
j));
"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;
"Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
err_div, div, t_div_from_op);
delta_diff2_val(
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)
"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[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
auto dm = simple_interface->
getDM();
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(
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(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
[](double, double, double) { return 1; })
);
};
auto solve_problem = [&] {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
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>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
auto base_mass = boost::make_shared<MatrixDouble>();
auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
base_mass, data_l2, base,
HCURL));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_grad));
pipeline_mng->getOpDomainRhsPipeline().push_back(
"FIELD1", ptr_divergence));
pipeline_mng->getOpDomainRhsPipeline().push_back(
ptr_hessian));
ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
CHKERR pipeline_mng->loopFiniteElements();
};
CHKERR assemble_matrices_and_vectors();
}
}
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.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(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)
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.
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....
default operator for TRI element
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.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
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.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
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