v0.13.0
Public Member Functions | Protected Attributes | List of all members
MyPostProc Struct Reference
Inheritance diagram for MyPostProc:
[legend]
Collaboration diagram for MyPostProc:
[legend]

Public Member Functions

MoFEMErrorCode generateReferenceElementMesh ()
 
MoFEMErrorCode setGaussPts (int order)
 

Protected Attributes

ublas::matrix< int > refEleMap
 
MatrixDouble shapeFunctions
 

Detailed Description

Examples
plot_base.cpp.

Definition at line 52 of file plot_base.cpp.

Member Function Documentation

◆ generateReferenceElementMesh()

MoFEMErrorCode MyPostProc::generateReferenceElementMesh ( )
Examples
plot_base.cpp.

Definition at line 354 of file plot_base.cpp.

354  {
356  moab::Core core_ref;
357  moab::Interface &moab_ref = core_ref;
358 
359  char ref_mesh_file_name[255];
360 
361  if (SPACE_DIM == 2)
362  strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
363  else if (SPACE_DIM == 3)
364  strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
365  else
366  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
367  "Dimension not implemented");
368 
369  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-ref_file", ref_mesh_file_name,
370  255, PETSC_NULL);
371  CHKERR moab_ref.load_file(ref_mesh_file_name, 0, "");
372 
373  // Get elements
374  Range elems;
375  CHKERR moab_ref.get_entities_by_dimension(0, SPACE_DIM, elems);
376 
377  // Add mid-nodes on edges
378  EntityHandle meshset;
379  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
380  CHKERR moab_ref.add_entities(meshset, elems);
381  CHKERR moab_ref.convert_entities(meshset, true, false, false);
382  CHKERR moab_ref.delete_entities(&meshset, 1);
383 
384  // Get nodes on the mesh
385  Range elem_nodes;
386  CHKERR moab_ref.get_connectivity(elems, elem_nodes, false);
387 
388  // Map node entity and Gauss pint number
389  std::map<EntityHandle, int> nodes_pts_map;
390 
391  // Set gauss points coordinates from the reference mesh
392  gaussPts.resize(SPACE_DIM + 1, elem_nodes.size(), false);
393  gaussPts.clear();
394  Range::iterator nit = elem_nodes.begin();
395  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
396  double coords[3];
397  CHKERR moab_ref.get_coords(&*nit, 1, coords);
398  for (auto d : {0, 1, 2})
399  gaussPts(d, gg) = coords[d];
400  nodes_pts_map[*nit] = gg;
401  }
402 
403  if (SPACE_DIM == 2) {
404  // Set size of adjacency matrix (note ho order nodes 3 nodes and 3 nodes on
405  // edges)
406  refEleMap.resize(elems.size(), 3 + 3);
407  } else if (SPACE_DIM == 3) {
408  refEleMap.resize(elems.size(), 4 + 6);
409  }
410 
411  // Set adjacency matrix
412  Range::iterator tit = elems.begin();
413  for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
414  const EntityHandle *conn;
415  int num_nodes;
416  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
417  for (int nn = 0; nn != num_nodes; ++nn) {
418  refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
419  }
420  }
421 
423 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
constexpr int SPACE_DIM
Definition: plot_base.cpp:44
ublas::matrix< int > refEleMap
Definition: plot_base.cpp:59

◆ setGaussPts()

MoFEMErrorCode MyPostProc::setGaussPts ( int  order)
Examples
plot_base.cpp.

Definition at line 425 of file plot_base.cpp.

425  {
427 
428  const int num_nodes = gaussPts.size2();
429 
430  // Calculate sheape functions
431 
432  switch (numeredEntFiniteElementPtr->getEntType()) {
433  case MBTRI:
434  shapeFunctions.resize(num_nodes, 3);
435  CHKERR Tools::shapeFunMBTRI(&*shapeFunctions.data().begin(),
436  &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
437  break;
438  case MBQUAD: {
439  shapeFunctions.resize(num_nodes, 4);
440  for (int gg = 0; gg != num_nodes; gg++) {
441  double ksi = gaussPts(0, gg);
442  double eta = gaussPts(1, gg);
443  shapeFunctions(gg, 0) = N_MBQUAD0(ksi, eta);
444  shapeFunctions(gg, 1) = N_MBQUAD1(ksi, eta);
445  shapeFunctions(gg, 2) = N_MBQUAD2(ksi, eta);
446  shapeFunctions(gg, 3) = N_MBQUAD3(ksi, eta);
447  }
448  } break;
449  case MBTET: {
450  shapeFunctions.resize(num_nodes, 8);
451  CHKERR Tools::shapeFunMBTET(&*shapeFunctions.data().begin(),
452  &gaussPts(0, 0), &gaussPts(1, 0),
453  &gaussPts(2, 0), num_nodes);
454  } break;
455  case MBHEX: {
456  shapeFunctions.resize(num_nodes, 8);
457  for (int gg = 0; gg != num_nodes; gg++) {
458  double ksi = gaussPts(0, gg);
459  double eta = gaussPts(1, gg);
460  double zeta = gaussPts(2, gg);
461  shapeFunctions(gg, 0) = N_MBHEX0(ksi, eta, zeta);
462  shapeFunctions(gg, 1) = N_MBHEX1(ksi, eta, zeta);
463  shapeFunctions(gg, 2) = N_MBHEX2(ksi, eta, zeta);
464  shapeFunctions(gg, 3) = N_MBHEX3(ksi, eta, zeta);
465  shapeFunctions(gg, 4) = N_MBHEX4(ksi, eta, zeta);
466  shapeFunctions(gg, 5) = N_MBHEX5(ksi, eta, zeta);
467  shapeFunctions(gg, 6) = N_MBHEX6(ksi, eta, zeta);
468  shapeFunctions(gg, 7) = N_MBHEX7(ksi, eta, zeta);
469  }
470  } break;
471  default:
472  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
473  "Not implemented element type");
474  }
475 
476 
477  // Create physical nodes
478 
479  // MoAB interface allowing for creating nodes and elements in the bulk
480  ReadUtilIface *iface;
481  CHKERR postProcMesh.query_interface(iface);
482 
483  std::vector<double *> arrays; /// pointers to memory allocated by MoAB for
484  /// storing X, Y, and Z coordinates
485  EntityHandle startv; // Starting handle for vertex
486  // Allocate memory for num_nodes, and return starting handle, and access to
487  // memort.
488  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
489 
490  mapGaussPts.resize(gaussPts.size2());
491  for (int gg = 0; gg != num_nodes; ++gg)
492  mapGaussPts[gg] = startv + gg;
493 
494  Tag th;
495  int def_in_the_loop = -1;
496  CHKERR postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th,
497  MB_TAG_CREAT | MB_TAG_SPARSE,
498  &def_in_the_loop);
499 
500  // Create physical elements
501 
502  const int num_el = refEleMap.size1();
503  const int num_nodes_on_ele = refEleMap.size2();
504 
505  EntityHandle starte; // Starting handle to first created element
506  EntityHandle *conn; // Access to MOAB memory with connectivity of elements
507 
508  // Create tris/tets in the bulk in MoAB database
509  if (SPACE_DIM == 2)
510  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
511  starte, conn);
512  else if (SPACE_DIM == 3)
513  CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
514  starte, conn);
515  else
516  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
517  "Dimension not implemented");
518 
519  // At this point elements (memory for elements) is allocated, at code bellow
520  // actual connectivity of elements is set.
521  for (unsigned int tt = 0; tt != refEleMap.size1(); ++tt) {
522  for (int nn = 0; nn != num_nodes_on_ele; ++nn)
523  conn[num_nodes_on_ele * tt + nn] = mapGaussPts[refEleMap(tt, nn)];
524  }
525 
526  // Finalise elements creation. At that point MOAB updates adjacency tables,
527  // and elements are ready to use.
528  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
529 
530  auto physical_elements = Range(starte, starte + num_el - 1);
531  CHKERR postProcMesh.tag_clear_data(th, physical_elements, &(nInTheLoop));
532 
533  EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
534  int fe_num_nodes;
535  {
536  const EntityHandle *conn;
537  mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes, true);
538  coords.resize(3 * fe_num_nodes, false);
539  CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
540  }
541 
542  // Set physical coordinates to physical nodes
545  &*shapeFunctions.data().begin());
546 
548  arrays[0], arrays[1], arrays[2]);
549  const double *t_coords_ele_x = &coords[0];
550  const double *t_coords_ele_y = &coords[1];
551  const double *t_coords_ele_z = &coords[2];
552  for (int gg = 0; gg != num_nodes; ++gg) {
554  t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
555  t_coords(i) = 0;
556  for (int nn = 0; nn != fe_num_nodes; ++nn) {
557  t_coords(i) += t_n * t_ele_coords(i);
558  for (auto ii : {0, 1, 2})
559  if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
560  t_coords(ii) = 0;
561  ++t_ele_coords;
562  ++t_n;
563  }
564  ++t_coords;
565  }
566 
568 }
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:70
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:88
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:84
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:86
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:85
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:81
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:87
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:83
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:67
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:82
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:69
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:68
constexpr double eta
FTensor::Index< 'i', SPACE_DIM > i
MatrixDouble shapeFunctions
Definition: plot_base.cpp:60

Member Data Documentation

◆ refEleMap

ublas::matrix<int> MyPostProc::refEleMap
protected

Definition at line 59 of file plot_base.cpp.

◆ shapeFunctions

MatrixDouble MyPostProc::shapeFunctions
protected

Definition at line 60 of file plot_base.cpp.


The documentation for this struct was generated from the following file: