Inheritance diagram for Op< OP >:

Collaboration diagram for Op< OP >:

Public Member Functions
	Op (moab::Interface &post_proc, MapCoords &map_coords, EntityHandle prism)

MoFEMErrorCode	doWork (int side, EntityType type, EntitiesFieldData::EntData &data)

Public Attributes
moab::Interface &	postProc

MapCoords &	mapCoords

EntityHandle	prism

std::vector< EntityHandle >	nodeHandles

Detailed Description

template<typename OP>
struct Op< OP >

Examples: build_large_problem.cpp, and prism_elements_from_surface.cpp.

Definition at line 78 of file prism_elements_from_surface.cpp.

Constructor & Destructor Documentation

◆ Op()

template<typename OP >

Op< OP >::Op	(	moab::Interface &	post_proc,
		MapCoords &	map_coords,
		EntityHandle	prism
	)

Examples: prism_elements_from_surface.cpp.

Definition at line 579 of file prism_elements_from_surface.cpp.

581 : OP("FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),

582 postProc(post_proc), mapCoords(map_coords), prism(prism) {}

OP

Op::mapCoords

MapCoords & mapCoords

Definition: prism_elements_from_surface.cpp:86

Op::postProc

moab::Interface & postProc

Definition: prism_elements_from_surface.cpp:85

Op::prism

EntityHandle prism

Definition: prism_elements_from_surface.cpp:87

Member Function Documentation

◆ doWork()

template<typename OP >

MoFEMErrorCode Op< OP >::doWork	(	int	side,
		EntityType	type,
		EntitiesFieldData::EntData &	data
	)

Examples: build_large_problem.cpp, and prism_elements_from_surface.cpp.

Definition at line 585 of file prism_elements_from_surface.cpp.

                                                              {
  constexpr double def_val[] = {0, 0, 0};
  MoFEMFunctionBegin;
  switch (type) {
  case MBVERTEX:
  case MBEDGE:
  case MBTRI:
  case MBQUAD:
    break;
  default:
    MoFEMFunctionReturnHot(0);
  }
 
  const int nb_dofs = data.getIndices().size();
  for (int dd = 0; dd != nb_dofs; ++dd)
    if (data.getFieldData()[dd] != data.getIndices()[dd]) {
      std::cerr << "Indices: " << data.getIndices() << std::endl;
      std::cerr << "Local indices: " << data.getLocalIndices() << std::endl;
      std::cerr << "Data: " << data.getFieldData() << std::endl;
      SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Indicices/data inconsistency %3.1f != %d",
               data.getFieldData()[dd], data.getIndices()[dd]);
    }
 
  const EntityHandle fe_ent = OP::getFEEntityHandle();
  const EntityHandle ent = OP::getSideEntity(side, type);
  int side_prism, sense, offset;
  if (type == MBVERTEX) {
    CHKERR postProc.side_number(prism, fe_ent, side_prism, sense, offset);
  } else
    CHKERR postProc.side_number(prism, ent, side_prism, sense, offset);
 
  if (type == MBVERTEX) {
    auto &coords_at_pts = OP::getCoordsAtGaussPts();
    const size_t nb_gauss_pts = coords_at_pts.size1();
 
    nodeHandles.reserve(nb_gauss_pts);
    nodeHandles.clear();
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
      auto t = boost::make_tuple(CoordsAndHandle::getArg(coords_at_pts(gg, 0)),
                                 CoordsAndHandle::getArg(coords_at_pts(gg, 1)),
                                 CoordsAndHandle::getArg(coords_at_pts(gg, 2)));
 
      auto it = mapCoords.find(t);
      if (it != mapCoords.end())
        nodeHandles.emplace_back(it->node);
      else
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vertex not found");
    }
  }
 
  auto to_str = [](auto i) { return boost::lexical_cast<std::string>(i); };
  std::string tag_name_base =
      "FEType" + to_str(OP::getNumeredEntFiniteElementPtr()->getEntType()) +
      "Type" + to_str(type) + "Side" + to_str(side_prism);
 
  std::string tag_prism_name_base =
      "PrismType" + to_str(type) + "Side" + to_str(side_prism);
 
  MatrixDouble trans_base = trans(data.getN());
  MatrixDouble prism_base(trans_base.size1(), trans_base.size2());
  if (trans_base.size2() != nodeHandles.size())
    SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "wrong size %d != %d",
             trans_base.size2(), nodeHandles.size());
 
  for (size_t rr = 0; rr != trans_base.size1(); ++rr) {
 
    std::string tag_name = tag_name_base + "Base";
    if (type == MBVERTEX) {
      EntityHandle node = data.getFieldDofs()[rr]->getEnt();
      int prism_mode_side;
      CHKERR postProc.side_number(prism, node, prism_mode_side, sense, offset);
      tag_name += to_str(prism_mode_side);
    } else {
      tag_name += to_str(rr);
    }
 
    std::cout << tag_name << endl;
 
    Tag th;
    CHKERR postProc.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, th,
                                   MB_TAG_CREAT | MB_TAG_DENSE, def_val);
    CHKERR postProc.tag_set_data(th, &nodeHandles[0], nodeHandles.size(),
                                 &trans_base(rr, 0));
 
    if (type != MBVERTEX) {
      auto tag_prism_name = tag_prism_name_base + "Base" + to_str(rr);
      Tag th_prism;
      CHKERR postProc.tag_get_handle(tag_prism_name.c_str(), th_prism);
      CHKERR postProc.tag_get_data(th_prism, &nodeHandles[0],
                                   nodeHandles.size(), &prism_base(rr, 0));
    }
  }
 
  auto sum_matrix = [](MatrixDouble &m) {
    double s = 0;
    for (unsigned int ii = 0; ii < m.size1(); ii++) {
      for (unsigned int jj = 0; jj < m.size2(); jj++) {
        s += std::abs(m(ii, jj));
      }
    }
    return s;
  };
 
  if (type != MBVERTEX) {
    prism_base -= trans_base;
    double sum = sum_matrix(prism_base);
    constexpr double eps = 1e-6;
 
    if (std::abs(sum) > eps)
      cout << "Inconsistent base " << tag_prism_name_base << " "
           << tag_name_base << " sum  " << sum << endl;
 
    if (std::abs(sum) > eps)
      SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Inconsistent base %s sum %6.4e", tag_prism_name_base.c_str(),
               sum);
  }
 
  MoFEMFunctionReturn(0);
}