262 {
264
267 pipeline_mng->getDomainLhsFE().reset();
268 pipeline_mng->getSkeletonRhsFE().reset();
269 pipeline_mng->getSkeletonLhsFE().reset();
270 pipeline_mng->getBoundaryRhsFE().reset();
271 pipeline_mng->getBoundaryLhsFE().reset();
272
273 auto rule = [](int, int, int p) -> int { return 2 * p; };
274 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
275 auto rule_2 = [
this](int, int, int) {
return 2 *
oRder; };
276 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
277
278 auto add_base_ops = [&](auto &pipeline) {
279 auto det_ptr = boost::make_shared<VectorDouble>();
280 auto jac_ptr = boost::make_shared<MatrixDouble>();
281 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
285 };
286
287 auto u_vals_ptr = boost::make_shared<VectorDouble>();
288 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
289
290 add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
293 pipeline_mng->getOpDomainRhsPipeline().push_back(
295
296 enum NORMS {
L2 = 0, SEMI_NORM,
H1, LAST_NORM };
297 std::array<int, LAST_NORM> error_indices;
298 for (
int i = 0;
i != LAST_NORM; ++
i)
299 error_indices[
i] =
i;
302
304 error_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side, EntityType
type,
308
310
311 if (const size_t nb_dofs = data.getIndices().size()) {
312
313 const int nb_integration_pts = o->getGaussPts().size2();
314 auto t_w = o->getFTensor0IntegrationWeight();
317 auto t_coords = o->getFTensor1CoordsAtGaussPts();
318
319 std::array<double, LAST_NORM> error;
320 std::fill(error.begin(), error.end(), 0);
321
322 for (int gg = 0; gg != nb_integration_pts; ++gg) {
323
324 const double alpha = t_w * o->getMeasure();
325 const double diff =
326 t_val -
u_exact(t_coords(0), t_coords(1), t_coords(2));
327
328 auto t_exact_grad =
u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
329
330 const double diff_grad =
331 (t_grad(
i) - t_exact_grad(
i)) * (t_grad(
i) - t_exact_grad(
i));
332
333 error[
L2] += alpha * pow(diff, 2);
334 error[SEMI_NORM] += alpha * diff_grad;
335
336 ++t_w;
337 ++t_val;
338 ++t_grad;
339 ++t_coords;
340 }
341
342 error[
H1] = error[
L2] + error[SEMI_NORM];
343
345 ADD_VALUES);
346 }
347
349 };
350
351 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
352 add_base_ops(side_fe_ptr->getOpPtrVector());
353 side_fe_ptr->getOpPtrVector().push_back(
355 std::array<VectorDouble, 2> side_vals;
356 std::array<double, 2> area_map;
357
359 side_vals_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
364
365 const auto nb_in_loop = o->getFEMethod()->nInTheLoop;
366 area_map[nb_in_loop] = o->getMeasure();
367 side_vals[nb_in_loop] = *u_vals_ptr;
368 if (!nb_in_loop) {
369 area_map[1] = 0;
370 side_vals[1].clear();
371 }
372
374 };
375 side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
376
377 auto do_work_rhs_error = [&](
DataOperator *op_ptr,
int side, EntityType
type,
381
382 CHKERR o->loopSideFaces(
"dFE", *side_fe_ptr);
383 const auto in_the_loop = side_fe_ptr->nInTheLoop;
384
385#ifndef NDEBUG
386 const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
388 << "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
389#endif
390
391 const double s = o->getMeasure() / (area_map[0] + area_map[1]);
393
394 std::array<double, LAST_NORM> error;
395 std::fill(error.begin(), error.end(), 0);
396
397 const int nb_integration_pts = o->getGaussPts().size2();
398
399 if (!in_the_loop) {
400 side_vals[1].resize(nb_integration_pts, false);
401 auto t_coords = o->getFTensor1CoordsAtGaussPts();
403 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
404 t_val_m =
u_exact(t_coords(0), t_coords(1), t_coords(2));
405 ++t_coords;
406 ++t_val_m;
407 }
408 }
409
412 auto t_w = o->getFTensor0IntegrationWeight();
413
414 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
415
416 const double alpha = o->getMeasure() * t_w;
417 const auto diff = t_val_p - t_val_m;
418 error[SEMI_NORM] += alpha * p * diff * diff;
419
420 ++t_w;
421 ++t_val_p;
422 ++t_val_m;
423 }
424
425
426 error[
H1] = error[SEMI_NORM];
428 ADD_VALUES);
429
431 };
432
434 skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
436 boundary_error_op->doWorkRhsHook = do_work_rhs_error;
437
438 pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
439 pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
440 pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
441
442 CHKERR pipeline_mng->loopFiniteElements();
443
444 CHKERR VecAssemblyBegin(l2_vec);
445 CHKERR VecAssemblyEnd(l2_vec);
446
448 const double *array;
449 CHKERR VecGetArrayRead(l2_vec, &array);
450 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm L2 %6.4e",
451 std::sqrt(array[
L2]));
452 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm Energetic %6.4e",
453 std::sqrt(array[SEMI_NORM]));
454 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm H1 %6.4e",
455 std::sqrt(array[
H1]));
456
458 constexpr double eps = 1e-12;
459 if (std::sqrt(array[
H1]) >
eps)
461 }
462
463 CHKERR VecRestoreArrayRead(l2_vec, &array);
468 }
469
470
471
473}
#define MOFEM_LOG_C(channel, severity, format,...)
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
BoundaryEle::UserDataOperator BoundaryEleOp
@ 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.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
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.
FTensor::Index< 'i', SPACE_DIM > i
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)
Get field gradients at integration pts for scalar field 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.