27static char help[] =
"...\n\n";
47 GAUSS>::OpSource<BASE_DIM, BASE_DIM>;
49constexpr double a0 = 0.0;
50constexpr double a1 = 2.0;
51constexpr double a2 = -15.0 *
a0;
52constexpr double a3 = -20.0 / 6 *
a1;
53constexpr double a4 = 15.0 *
a0;
62 6 *
a6 * std::pow(x, 5) * std::pow(y, 0) +
63 5 *
a5 * std::pow(x, 4) * std::pow(y, 1) +
64 4 *
a4 * std::pow(x, 3) * std::pow(y, 2) +
65 3 *
a3 * std::pow(x, 2) * std::pow(y, 3) +
66 2 *
a2 * std::pow(x, 1) * std::pow(y, 4) +
67 1 *
a1 * std::pow(x, 0) * std::pow(y, 5),
69 1 *
a5 * std::pow(x, 5) * std::pow(y, 0) +
70 2 *
a4 * std::pow(x, 4) * std::pow(y, 1) +
71 3 *
a3 * std::pow(x, 3) * std::pow(y, 2) +
72 4 *
a2 * std::pow(x, 2) * std::pow(y, 3) +
73 5 *
a1 * std::pow(x, 1) * std::pow(y, 4) +
74 6 *
a0 * std::pow(x, 0) * std::pow(y, 5),
84 30 *
a6 * pow(x, 4) * pow(y, 0) + 20 *
a5 * pow(x, 3) * pow(y, 1) +
85 12 *
a4 * pow(x, 2) * pow(y, 2) + 6 *
a3 * pow(x, 1) * pow(y, 3) +
86 2 *
a2 * pow(x, 0) * pow(y, 4),
88 5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
89 9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
90 5 *
a1 * pow(x, 0) * pow(y, 4),
92 5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
93 9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
94 5 *
a1 * pow(x, 0) * pow(y, 4),
96 2 *
a4 * pow(x, 4) * pow(y, 0) + 6 *
a3 * pow(x, 3) * pow(y, 1) +
97 12 *
a2 * pow(x, 2) * pow(y, 2) + 20 *
a1 * pow(x, 1) * pow(y, 3) +
98 30 *
a0 * pow(x, 0) * pow(y, 4),
108 30 * 4 *
a6 * pow(x, 3) * pow(y, 0) +
109 20 * 3 *
a5 * pow(x, 2) * pow(y, 1) +
110 12 * 2 *
a4 * pow(x, 1) * pow(y, 2) +
111 6 * 1 *
a3 * pow(x, 0) * pow(y, 3),
115 20 * 1 *
a5 * pow(x, 3) * pow(y, 0) +
116 12 * 2 *
a4 * pow(x, 2) * pow(y, 2) +
117 6 * 3 *
a3 * pow(x, 1) * pow(y, 2) +
118 2 * 4 *
a2 * pow(x, 0) * pow(y, 3),
122 5 * 4 *
a5 * pow(x, 3) * pow(y, 0) +
123 8 * 3 *
a4 * pow(x, 2) * pow(y, 1) +
124 9 * 2 *
a3 * pow(x, 1) * pow(y, 2) +
125 8 * 1 *
a2 * pow(x, 0) * pow(y, 3),
129 8 * 1 *
a4 * pow(x, 3) * pow(y, 0) +
130 9 * 2 *
a3 * pow(x, 2) * pow(y, 1) +
131 8 * 3 *
a2 * pow(x, 1) * pow(y, 2) +
132 5 * 4 *
a1 * pow(x, 0) * pow(y, 3),
136 5 * 4 *
a5 * pow(x, 3) * pow(y, 0) +
137 8 * 3 *
a4 * pow(x, 2) * pow(y, 1) +
138 9 * 2 *
a3 * pow(x, 1) * pow(y, 2) +
139 8 * 1 *
a2 * pow(x, 0) * pow(y, 3),
143 8 * 1 *
a4 * pow(x, 3) * pow(y, 0) +
144 9 * 2 *
a3 * pow(x, 2) * pow(y, 1) +
145 8 * 3 *
a2 * pow(x, 1) * pow(y, 2) +
146 5 * 4 *
a1 * pow(x, 0) * pow(y, 3),
150 2 * 4 *
a4 * pow(x, 3) * pow(y, 0) +
151 6 * 3 *
a3 * pow(x, 2) * pow(y, 1) +
152 12 * 2 *
a2 * pow(x, 1) * pow(y, 2) +
153 20 * 1 *
a1 * pow(x, 0) * pow(y, 3),
157 6 * 1 *
a3 * pow(x, 3) * pow(y, 0) +
158 12 * 2 *
a2 * pow(x, 2) * pow(y, 1) +
159 20 * 3 *
a1 * pow(x, 1) * pow(y, 2) +
160 30 * 4 *
a0 * pow(x, 0) * pow(y, 3),
180 boost::shared_ptr<VectorDouble> ptr_div,
181 boost::shared_ptr<MatrixDouble> ptr_grad,
182 boost::shared_ptr<MatrixDouble> ptr_hess)
193 const double eps = 1e-6;
194 if (
type == MBEDGE && side == 0) {
195 const int nb_gauss_pts = data.
getN().size1();
197 auto t_vals_from_op = getFTensor1FromMat<3>(*
ptrVals);
199 auto t_grad_from_op = getFTensor2FromMat<3, 2>(*
ptrGrad);
200 auto t_hess_from_op = getFTensor3FromMat<3, 2, 2>(*
ptrHess);
202 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
209 double err_val = sqrt(delta_val(
i) * delta_val(
i));
212 "Wrong value %4.3e", err_val);
215 delta_diff_val(
i,
j) =
217 double err_diff_val = sqrt(delta_diff_val(
i,
j) * delta_diff_val(
i,
j));
218 if (err_diff_val >
eps)
220 "Wrong derivative of value %4.3e", err_diff_val);
222 double div = t_grad_from_op(0, 0) + t_grad_from_op(1, 1);
223 double err_div = div - t_div_from_op;
226 "Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
227 err_div, div, t_div_from_op);
230 delta_diff2_val(
i,
j,
k) =
232 double hess_diff_error =
233 sqrt(delta_diff2_val(
i,
j,
k) * delta_diff2_val(
i,
j,
k));
234 if (hess_diff_error >
eps)
236 "Wrong hessian from operator %4.3e", hess_diff_error);
248int main(
int argc,
char *argv[]) {
254 DMType dm_name =
"DMMOFEM";
270 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
271 const char *list_bases[] = {
"ainsworth",
"demkowicz"};
273 PetscInt choice_base_value = AINSWORTH;
275 LASBASETOP, &choice_base_value, &flg);
277 if (flg != PETSC_TRUE)
280 if (choice_base_value == AINSWORTH)
282 else if (choice_base_value == DEMKOWICZ)
286 constexpr int order = 5;
289 auto dm = simple_interface->
getDM();
293 auto assemble_matrices_and_vectors = [&]() {
295 auto jac_ptr = boost::make_shared<MatrixDouble>();
296 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
297 auto det_ptr = boost::make_shared<VectorDouble>();
299 pipeline_mng->getOpDomainRhsPipeline().push_back(
301 pipeline_mng->getOpDomainRhsPipeline().push_back(
303 pipeline_mng->getOpDomainRhsPipeline().push_back(
305 pipeline_mng->getOpDomainRhsPipeline().push_back(
307 pipeline_mng->getOpDomainRhsPipeline().push_back(
310 pipeline_mng->getOpDomainLhsPipeline().push_back(
312 pipeline_mng->getOpDomainLhsPipeline().push_back(
314 pipeline_mng->getOpDomainLhsPipeline().push_back(
316 pipeline_mng->getOpDomainLhsPipeline().push_back(
318 pipeline_mng->getOpDomainLhsPipeline().push_back(
321 [](
double,
double,
double) {
return 1; })
332 auto solve_problem = [&] {
334 auto solver = pipeline_mng->createKSP();
335 CHKERR KSPSetFromOptions(solver);
342 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
343 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
348 auto check_solution = [&] {
351 pipeline_mng->getOpDomainLhsPipeline().clear();
352 pipeline_mng->getOpDomainRhsPipeline().clear();
354 auto ptr_values = boost::make_shared<MatrixDouble>();
355 auto ptr_divergence = boost::make_shared<VectorDouble>();
356 auto ptr_grad = boost::make_shared<MatrixDouble>();
357 auto ptr_hessian = boost::make_shared<MatrixDouble>();
359 auto jac_ptr = boost::make_shared<MatrixDouble>();
360 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
361 auto det_ptr = boost::make_shared<VectorDouble>();
364 pipeline_mng->getOpDomainRhsPipeline().push_back(
366 pipeline_mng->getOpDomainRhsPipeline().push_back(
368 pipeline_mng->getOpDomainRhsPipeline().push_back(
370 pipeline_mng->getOpDomainRhsPipeline().push_back(
372 pipeline_mng->getOpDomainRhsPipeline().push_back(
376 auto base_mass = boost::make_shared<MatrixDouble>();
377 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
379 pipeline_mng->getOpDomainRhsPipeline().push_back(
381 pipeline_mng->getOpDomainRhsPipeline().push_back(
384 pipeline_mng->getOpDomainRhsPipeline().push_back(
386 base_mass, data_l2, base,
HCURL));
389 pipeline_mng->getOpDomainRhsPipeline().push_back(
392 pipeline_mng->getOpDomainRhsPipeline().push_back(
396 pipeline_mng->getOpDomainRhsPipeline().push_back(
398 "FIELD1", ptr_divergence));
399 pipeline_mng->getOpDomainRhsPipeline().push_back(
404 ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
406 CHKERR pipeline_mng->loopFiniteElements();
411 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 smartCreateDMVector(DM dm)
Get smart vector from DM.
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.
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 > smartVectorDuplicate(SmartPetscObj< 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)
DeprecatedCoreInterface Interface
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)
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
friend class UserDataOperator
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 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 loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
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