15static char help[] =
"...\n\n";
35 GAUSS>::OpSource<BASE_DIM, BASE_DIM>;
37constexpr double a0 = 0.0;
38constexpr double a1 = 2.0;
39constexpr double a2 = -15.0 *
a0;
40constexpr double a3 = -20.0 / 6 *
a1;
41constexpr double a4 = 15.0 *
a0;
50 6 *
a6 * std::pow(x, 5) * std::pow(y, 0) +
51 5 *
a5 * std::pow(x, 4) * std::pow(y, 1) +
52 4 *
a4 * std::pow(x, 3) * std::pow(y, 2) +
53 3 *
a3 * std::pow(x, 2) * std::pow(y, 3) +
54 2 *
a2 * std::pow(x, 1) * std::pow(y, 4) +
55 1 *
a1 * std::pow(x, 0) * std::pow(y, 5),
57 1 *
a5 * std::pow(x, 5) * std::pow(y, 0) +
58 2 *
a4 * std::pow(x, 4) * std::pow(y, 1) +
59 3 *
a3 * std::pow(x, 3) * std::pow(y, 2) +
60 4 *
a2 * std::pow(x, 2) * std::pow(y, 3) +
61 5 *
a1 * std::pow(x, 1) * std::pow(y, 4) +
62 6 *
a0 * std::pow(x, 0) * std::pow(y, 5),
72 30 *
a6 * pow(x, 4) * pow(y, 0) + 20 *
a5 * pow(x, 3) * pow(y, 1) +
73 12 *
a4 * pow(x, 2) * pow(y, 2) + 6 *
a3 * pow(x, 1) * pow(y, 3) +
74 2 *
a2 * pow(x, 0) * pow(y, 4),
76 5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
77 9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
78 5 *
a1 * pow(x, 0) * pow(y, 4),
80 5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
81 9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
82 5 *
a1 * pow(x, 0) * pow(y, 4),
84 2 *
a4 * pow(x, 4) * pow(y, 0) + 6 *
a3 * pow(x, 3) * pow(y, 1) +
85 12 *
a2 * pow(x, 2) * pow(y, 2) + 20 *
a1 * pow(x, 1) * pow(y, 3) +
86 30 *
a0 * pow(x, 0) * pow(y, 4),
96 30 * 4 *
a6 * pow(x, 3) * pow(y, 0) +
97 20 * 3 *
a5 * pow(x, 2) * pow(y, 1) +
98 12 * 2 *
a4 * pow(x, 1) * pow(y, 2) +
99 6 * 1 *
a3 * pow(x, 0) * pow(y, 3),
103 20 * 1 *
a5 * pow(x, 3) * pow(y, 0) +
104 12 * 2 *
a4 * pow(x, 2) * pow(y, 2) +
105 6 * 3 *
a3 * pow(x, 1) * pow(y, 2) +
106 2 * 4 *
a2 * pow(x, 0) * pow(y, 3),
110 5 * 4 *
a5 * pow(x, 3) * pow(y, 0) +
111 8 * 3 *
a4 * pow(x, 2) * pow(y, 1) +
112 9 * 2 *
a3 * pow(x, 1) * pow(y, 2) +
113 8 * 1 *
a2 * pow(x, 0) * pow(y, 3),
117 8 * 1 *
a4 * pow(x, 3) * pow(y, 0) +
118 9 * 2 *
a3 * pow(x, 2) * pow(y, 1) +
119 8 * 3 *
a2 * pow(x, 1) * pow(y, 2) +
120 5 * 4 *
a1 * pow(x, 0) * pow(y, 3),
124 5 * 4 *
a5 * pow(x, 3) * pow(y, 0) +
125 8 * 3 *
a4 * pow(x, 2) * pow(y, 1) +
126 9 * 2 *
a3 * pow(x, 1) * pow(y, 2) +
127 8 * 1 *
a2 * pow(x, 0) * pow(y, 3),
131 8 * 1 *
a4 * pow(x, 3) * pow(y, 0) +
132 9 * 2 *
a3 * pow(x, 2) * pow(y, 1) +
133 8 * 3 *
a2 * pow(x, 1) * pow(y, 2) +
134 5 * 4 *
a1 * pow(x, 0) * pow(y, 3),
138 2 * 4 *
a4 * pow(x, 3) * pow(y, 0) +
139 6 * 3 *
a3 * pow(x, 2) * pow(y, 1) +
140 12 * 2 *
a2 * pow(x, 1) * pow(y, 2) +
141 20 * 1 *
a1 * pow(x, 0) * pow(y, 3),
145 6 * 1 *
a3 * pow(x, 3) * pow(y, 0) +
146 12 * 2 *
a2 * pow(x, 2) * pow(y, 1) +
147 20 * 3 *
a1 * pow(x, 1) * pow(y, 2) +
148 30 * 4 *
a0 * pow(x, 0) * pow(y, 3),
168 boost::shared_ptr<VectorDouble> ptr_div,
169 boost::shared_ptr<MatrixDouble> ptr_grad,
170 boost::shared_ptr<MatrixDouble> ptr_hess)
181 const double eps = 1e-6;
182 if (type == MBEDGE && side == 0) {
183 const int nb_gauss_pts = data.
getN().size1();
185 auto t_vals_from_op = getFTensor1FromMat<3>(*
ptrVals);
187 auto t_grad_from_op = getFTensor2FromMat<3, 2>(*
ptrGrad);
188 auto t_hess_from_op = getFTensor3FromMat<3, 2, 2>(*
ptrHess);
190 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
196 delta_val(
i) = t_vals_from_op(
i) - ApproxFunctions::fUn(x, y, 0)(
i);
197 double err_val = sqrt(delta_val(
i) * delta_val(
i));
200 "Wrong value %4.3e", err_val);
203 delta_diff_val(
i,
j) =
204 t_grad_from_op(
i,
j) - ApproxFunctions::diffFun(x, y)(
i,
j);
205 double err_diff_val = sqrt(delta_diff_val(
i,
j) * delta_diff_val(
i,
j));
206 if (err_diff_val >
eps)
208 "Wrong derivative of value %4.3e", err_diff_val);
210 double div = t_grad_from_op(0, 0) + t_grad_from_op(1, 1);
211 double err_div = div - t_div_from_op;
214 "Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
215 err_div, div, t_div_from_op);
218 delta_diff2_val(
i,
j,
k) =
219 t_hess_from_op(
i,
j,
k) - ApproxFunctions::diff2Fun(x, y)(
i,
j,
k);
220 double hess_diff_error =
221 sqrt(delta_diff2_val(
i,
j,
k) * delta_diff2_val(
i,
j,
k));
222 if (hess_diff_error >
eps)
224 "Wrong hessian from operator %4.3e", hess_diff_error);
236int main(
int argc,
char *argv[]) {
242 DMType dm_name =
"DMMOFEM";
245 moab::Core mb_instance;
246 moab::Interface &moab = mb_instance;
258 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
259 const char *list_bases[] = {
"ainsworth",
"demkowicz"};
261 PetscInt choice_base_value = AINSWORTH;
263 LASBASETOP, &choice_base_value, &flg);
265 if (flg != PETSC_TRUE)
268 if (choice_base_value == AINSWORTH)
270 else if (choice_base_value == DEMKOWICZ)
274 constexpr int order = 5;
277 auto dm = simple_interface->
getDM();
281 auto assemble_matrices_and_vectors = [&]() {
283 auto jac_ptr = boost::make_shared<MatrixDouble>();
284 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
285 auto det_ptr = boost::make_shared<VectorDouble>();
287 pipeline_mng->getOpDomainRhsPipeline().push_back(
289 pipeline_mng->getOpDomainRhsPipeline().push_back(
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
293 pipeline_mng->getOpDomainRhsPipeline().push_back(
295 pipeline_mng->getOpDomainRhsPipeline().push_back(
298 pipeline_mng->getOpDomainLhsPipeline().push_back(
300 pipeline_mng->getOpDomainLhsPipeline().push_back(
302 pipeline_mng->getOpDomainLhsPipeline().push_back(
304 pipeline_mng->getOpDomainLhsPipeline().push_back(
306 pipeline_mng->getOpDomainLhsPipeline().push_back(
309 [](
double,
double,
double) {
return 1; })
320 auto solve_problem = [&] {
322 auto solver = pipeline_mng->createKSP();
323 CHKERR KSPSetFromOptions(solver);
330 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
331 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
336 auto check_solution = [&] {
339 pipeline_mng->getOpDomainLhsPipeline().clear();
340 pipeline_mng->getOpDomainRhsPipeline().clear();
342 auto ptr_values = boost::make_shared<MatrixDouble>();
343 auto ptr_divergence = boost::make_shared<VectorDouble>();
344 auto ptr_grad = boost::make_shared<MatrixDouble>();
345 auto ptr_hessian = boost::make_shared<MatrixDouble>();
347 auto jac_ptr = boost::make_shared<MatrixDouble>();
348 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
349 auto det_ptr = boost::make_shared<VectorDouble>();
352 pipeline_mng->getOpDomainRhsPipeline().push_back(
354 pipeline_mng->getOpDomainRhsPipeline().push_back(
356 pipeline_mng->getOpDomainRhsPipeline().push_back(
358 pipeline_mng->getOpDomainRhsPipeline().push_back(
360 pipeline_mng->getOpDomainRhsPipeline().push_back(
364 auto base_mass = boost::make_shared<MatrixDouble>();
365 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
367 pipeline_mng->getOpDomainRhsPipeline().push_back(
369 pipeline_mng->getOpDomainRhsPipeline().push_back(
372 pipeline_mng->getOpDomainRhsPipeline().push_back(
374 base_mass, data_l2, base,
HCURL));
377 pipeline_mng->getOpDomainRhsPipeline().push_back(
380 pipeline_mng->getOpDomainRhsPipeline().push_back(
384 pipeline_mng->getOpDomainRhsPipeline().push_back(
386 "FIELD1", ptr_divergence));
387 pipeline_mng->getOpDomainRhsPipeline().push_back(
392 ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
394 CHKERR pipeline_mng->loopFiniteElements();
399 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.
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
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.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data 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.
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