261 {
263
266 pipeline_mng->getDomainLhsFE().reset();
267 pipeline_mng->getSkeletonRhsFE().reset();
268 pipeline_mng->getSkeletonLhsFE().reset();
269 pipeline_mng->getBoundaryRhsFE().reset();
270 pipeline_mng->getBoundaryLhsFE().reset();
271
272 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
273 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
275 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
276
277 auto add_base_ops = [&](auto &pipeline) {
278 auto det_ptr = boost::make_shared<VectorDouble>();
279 auto jac_ptr = boost::make_shared<MatrixDouble>();
280 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
284 };
285
286 auto u_vals_ptr = boost::make_shared<VectorDouble>();
287 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
288
289 add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
290 pipeline_mng->getOpDomainRhsPipeline().push_back(
292 pipeline_mng->getOpDomainRhsPipeline().push_back(
294
295 enum NORMS {
L2 = 0, SEMI_NORM,
H1, LAST_NORM };
296 std::array<int, LAST_NORM> error_indices;
297 for (
int i = 0;
i != LAST_NORM; ++
i)
298 error_indices[
i] =
i;
301
307
309
310 if (const size_t nb_dofs = data.getIndices().size()) {
311
312 const int nb_integration_pts =
o->getGaussPts().size2();
313 auto t_w =
o->getFTensor0IntegrationWeight();
315 auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
316 auto t_coords =
o->getFTensor1CoordsAtGaussPts();
317
318 std::array<double, LAST_NORM> error;
319 std::fill(error.begin(), error.end(), 0);
320
321 for (int gg = 0; gg != nb_integration_pts; ++gg) {
322
323 const double alpha = t_w *
o->getMeasure();
324 const double diff =
325 t_val -
u_exact(t_coords(0), t_coords(1), t_coords(2));
326
327 auto t_exact_grad =
u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
328
329 const double diff_grad =
330 (t_grad(
i) - t_exact_grad(
i)) * (t_grad(
i) - t_exact_grad(
i));
331
332 error[
L2] +=
alpha * pow(diff, 2);
333 error[SEMI_NORM] +=
alpha * diff_grad;
334
335 ++t_w;
336 ++t_val;
337 ++t_grad;
338 ++t_coords;
339 }
340
341 error[
H1] = error[
L2] + error[SEMI_NORM];
342
344 ADD_VALUES);
345 }
346
348 };
349
350 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
351 add_base_ops(side_fe_ptr->getOpPtrVector());
352 side_fe_ptr->getOpPtrVector().push_back(
354 std::array<VectorDouble, 2> side_vals;
355 std::array<double, 2> area_map;
356
358 side_vals_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
363
364 const auto nb_in_loop =
o->getFEMethod()->nInTheLoop;
365 area_map[nb_in_loop] =
o->getMeasure();
366 side_vals[nb_in_loop] = *u_vals_ptr;
367 if (!nb_in_loop) {
368 area_map[1] = 0;
369 side_vals[1].clear();
370 }
371
373 };
374 side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
375
380
381 CHKERR o->loopSideFaces(
"dFE", *side_fe_ptr);
382 const auto in_the_loop = side_fe_ptr->nInTheLoop;
383
384#ifndef NDEBUG
385 const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
387 << "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
388#endif
389
390 const double s =
o->getMeasure() / (area_map[0] + area_map[1]);
392
393 constexpr std::array<int, 2> sign_array{1, -1};
394
395 std::array<double, LAST_NORM> error;
396 std::fill(error.begin(), error.end(), 0);
397
398 const int nb_integration_pts =
o->getGaussPts().size2();
399
400 if (!in_the_loop) {
401 side_vals[1].resize(nb_integration_pts, false);
402 auto t_coords =
o->getFTensor1CoordsAtGaussPts();
404 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
405 t_val_m =
u_exact(t_coords(0), t_coords(1), t_coords(2));
406 ++t_coords;
407 ++t_val_m;
408 }
409 }
410
413 auto t_w =
o->getFTensor0IntegrationWeight();
414
415 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
416
417 const double alpha =
o->getMeasure() * t_w;
418 const auto diff = t_val_p - t_val_m;
419 error[SEMI_NORM] +=
alpha *
p * diff * diff;
420
421 ++t_w;
422 ++t_val_p;
423 ++t_val_m;
424 }
425
426
427 error[
H1] = error[SEMI_NORM];
429 ADD_VALUES);
430
432 };
433
435 skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
437 boundary_error_op->doWorkRhsHook = do_work_rhs_error;
438
439 pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
440 pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
441 pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
442
443 CHKERR pipeline_mng->loopFiniteElements();
444
445 CHKERR VecAssemblyBegin(l2_vec);
446 CHKERR VecAssemblyEnd(l2_vec);
447
449 const double *array;
450 CHKERR VecGetArrayRead(l2_vec, &array);
451 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm L2 %6.4e",
452 std::sqrt(array[
L2]));
453 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm Energetic %6.4e",
454 std::sqrt(array[SEMI_NORM]));
455 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm H1 %6.4e",
456 std::sqrt(array[
H1]));
457
459 constexpr double eps = 1e-12;
460 if (std::sqrt(array[
H1]) >
eps)
462 }
463
464 CHKERR VecRestoreArrayRead(l2_vec, &array);
469 }
470
471
472
474}
#define MOFEM_LOG_C(channel, severity, format,...)
@ L2
field with C-1 continuity
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
DomainEle::UserDataOperator DomainEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
base operator to do operations at Gauss Pt. level
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
keeps basic data about problem
DofIdx getNbDofsRow() const
MoFEMErrorCode getDM(DM *dm)
Get DM.