Projection10NodeCoordsOnField.cpp

Go to the documentation of this file.

/** \file Projection10NodeCoordsOnField.cpp
 * \brief Project displacements/coordinates from 10 node tetrahedra on
 * hierarchical approximation base.
 */
 
 
 
namespace MoFEM {
 
Projection10NodeCoordsOnField::Projection10NodeCoordsOnField(
    Interface &m_field, std::string field_name, int verb)
    : mField(m_field), fieldName(field_name), vErbose(verb) {}
 
MoFEMErrorCode Projection10NodeCoordsOnField::preProcess() {
  MoFEMFunctionBeginHot;
  auto field_ptr = mField.get_field_structure(fieldName);
  if (field_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
    MOFEM_TAG_AND_LOG("WORLD", Sev::warning, "Projection10NodeCoordsOnField")
        << "Only working well for first order AINSWORTH_BERNSTEIN_BEZIER_BASE!";
    MOFEM_LOG_CHANNEL("WORLD");
  }
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode Projection10NodeCoordsOnField::operator()() {
  MoFEMFunctionBeginHot;
 
  if (dofPtr == NULL) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
  }
  if (dofPtr->getName() != fieldName)
    MoFEMFunctionReturnHot(0);
  if (dofPtr->getEntType() == MBVERTEX) {
    EntityHandle node = dofPtr->getEnt();
    coords.resize(3);
    CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
    dofPtr->getFieldData() = coords[dofPtr->getDofCoeffIdx()];
    if (vErbose > 0) {
      MOFEM_TAG_AND_LOG_C("SELF", Sev::noisy, "Projection10NodeCoordsOnField",
                          "val = %6.7e\n", dofPtr->getFieldData());
      MOFEM_LOG_CHANNEL("SELF");
    }
 
    MoFEMFunctionReturnHot(0);
  }
  if (dofPtr->getEntType() != MBEDGE) {
    MoFEMFunctionReturnHot(0);
  }
  if (dofPtr->getEntDofIdx() != dofPtr->getDofCoeffIdx()) {
    MoFEMFunctionReturnHot(0);
  }
  EntityHandle edge = dofPtr->getEnt();
  if (type_from_handle(edge) != MBEDGE) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only elements which are type of MBEDGE");
  }
  // coords
  int num_nodes;
  const EntityHandle *conn;
  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
  if ((num_nodes != 2) && (num_nodes != 3)) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only 4 node and 10 node tets");
  }
  coords.resize(num_nodes * 3);
  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
  aveMidCoord.resize(3);
  midNodeCoord.resize(3);
  for (int dd = 0; dd < 3; dd++) {
    aveMidCoord[dd] = (coords[0 * 3 + dd] + coords[1 * 3 + dd]) * 0.5;
    if (num_nodes == 3) {
      midNodeCoord[dd] = coords[2 * 3 + dd];
    } else {
      midNodeCoord[dd] = aveMidCoord[dd];
    }
  }
 
  const auto base = dofPtr->getApproxBase();
  double edge_shape_function_val;
  switch (base) {
  case AINSWORTH_LEGENDRE_BASE:
    edge_shape_function_val = 0.25;
    break;
  case AINSWORTH_LOBATTO_BASE:
    edge_shape_function_val = 0.25 * LOBATTO_PHI0(0);
    break;
  case DEMKOWICZ_JACOBI_BASE: {
    double L[3];
    CHKERR Legendre_polynomials(2, 0, NULL, L, NULL, 1);
    edge_shape_function_val = 0.125 * LOBATTO_PHI0(0);
  }; break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
  }
 
  diffNodeCoord.resize(3);
  ublas::noalias(diffNodeCoord) = midNodeCoord - aveMidCoord;
  dOf.resize(3);
  ublas::noalias(dOf) = diffNodeCoord / edge_shape_function_val;
  if (dofPtr->getNbOfCoeffs() > 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only fields which are rank 3 or lower");
  }
  dofPtr->getFieldData() = dOf[dofPtr->getDofCoeffIdx()];
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode Projection10NodeCoordsOnField::postProcess() {
  MoFEMFunctionBeginHot;
  MoFEMFunctionReturnHot(0);
}
 
// remove triangles which have at least free_edge_threshold edges which are on the outer boundary of the set
static MoFEMErrorCode filterOuterTrisByEdges(moab::Interface &moab,
                                             Range &tris,
                                             const int free_edge_threshold = 1) {
  MoFEMFunctionBeginHot;
 
  if (tris.empty()) {
    MoFEMFunctionReturnHot(0);
  }
 
  Range edges;
  CHKERR moab.get_adjacencies(tris, 1, true, edges, moab::Interface::UNION);
 
  std::unordered_map<EntityHandle, int> edge_adj_count;
  edge_adj_count.reserve(edges.size());
 
  for (auto edge : edges) {
    Range adj_tris;
    CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_tris,
                                moab::Interface::UNION);
    adj_tris = intersect(adj_tris, tris);
    edge_adj_count[edge] = static_cast<int>(adj_tris.size());
  }
 
  Range remove_tris;
  for (auto tri : tris) {
    Range tri_edges;
    CHKERR moab.get_adjacencies(&tri, 1, 1, false, tri_edges,
                                moab::Interface::UNION);
    int free_edges = 0;
    for (auto edge : tri_edges) {
      auto it = edge_adj_count.find(edge);
      if (it != edge_adj_count.end() && it->second == 1) {
        if (++free_edges >= free_edge_threshold) {
          remove_tris.insert(tri);
          break;
        }
      }
    }
  }
 
  if (!remove_tris.empty()) {
    tris = subtract(tris, remove_tris);
  }
 
  MoFEMFunctionReturnHot(0);
}
 
// remove triangles which have less than 2 adjacent triangles in the set, i.e. are on the outer boundary of the set
static MoFEMErrorCode filterOuterTrisByTris(moab::Interface &moab,
                                                Range &tris) {
  MoFEMFunctionBeginHot;
 
  if (tris.empty()) {
    MoFEMFunctionReturnHot(0);
  }
 
  Range remove_tris;
  for (auto tri : tris) {
 
    Range adj_edges;
    CHKERR moab.get_adjacencies(&tri, 1, 1, false, adj_edges,
                                moab::Interface::UNION);
    Range adj_tris;
    CHKERR moab.get_adjacencies(adj_edges, 2, false, adj_tris,
                                moab::Interface::UNION);
    adj_tris = intersect(adj_tris, tris);
    adj_tris.erase(tri);
    if (adj_tris.size() <= 1) {
      remove_tris.insert(tri);
    }
  }
 
  if (!remove_tris.empty()) {
    tris = subtract(tris, remove_tris);
  }
 
  MoFEMFunctionReturnHot(0);
}
 
Projection10NodeCoordsFromFineSurfaceOnField::
    Projection10NodeCoordsFromFineSurfaceOnField(Interface &m_field,
                                                 moab::Interface &fine_moab,
                                                 const Range &fine_tris,
                                                 std::string field_name,
                                                 double planar_angle_deg,
                                                 double max_disp_factor,
                                                 bool debug_surfaces,
                                                 int verb)
    : mField(m_field), fineMoab(fine_moab), fineTris(fine_tris),
      fieldName(field_name), vErbose(verb), fineTree(&fine_moab), fineRoot(0),
      fineFeatureTree(&fine_moab), fineFeatureRoot(0),
      planarAngleCos(std::cos(planar_angle_deg * 0.017453292519943295)),
      maxDispFactor(max_disp_factor), debugSurfaces(debug_surfaces) {}
 
MoFEMErrorCode Projection10NodeCoordsFromFineSurfaceOnField::preProcess() {
  MoFEMFunctionBeginHot;
 
  // Build coarse skin triangles from volume mesh.
  Range coarse_vols;
  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, coarse_vols);
  if (coarse_vols.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no 3d entities found in coarse mesh");
  }
 
  Skinner coarse_skinner(&mField.get_moab());
  coarseSkinTris.clear();
  CHKERR coarse_skinner.find_skin(0, coarse_vols, false, coarseSkinTris);
  if (coarseSkinTris.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no skin triangles found in coarse mesh");
  }
 
  coarseSkinEdges.clear();
  CHKERR mField.get_moab().get_adjacencies(coarseSkinTris, 1, false,
                                          coarseSkinEdges,
                                          moab::Interface::UNION);
 
  if (fineTris.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no surface triangles found in fine mesh");
  }
 
  auto tri_normal = [&](const EntityHandle tri,
                        FTensor::Tensor1<double, 3> &normal) {
    MoFEMFunctionBeginHot;
    const EntityHandle *conn;
    int num_nodes;
    double coords[9];
    CHKERR fineMoab.get_connectivity(tri, conn, num_nodes, true);
    CHKERR fineMoab.get_coords(conn, num_nodes, coords);
 
    CHKERR mField.getInterface<Tools>()->getTriNormal(coords, &normal(0));
    normal.normalize();
    MoFEMFunctionReturnHot(0);
  };
 
  const double fine_planar_cos =
      std::cos(std::acos(planarAngleCos) * 1e-6);
 
  Range fine_edges;
  CHKERR fineMoab.get_adjacencies(fineTris, 1, true, fine_edges,
                                  moab::Interface::UNION);
 
  if (debugSurfaces) {
    CHKERR CutMeshInterface::SaveData(fineMoab)("fine_edge_features.vtk",
                                                fine_edges);
  }
  fineFeatureTris.clear();
  for (auto edge : fine_edges) {
    Range adj_tris;
    CHKERR fineMoab.get_adjacencies(&edge, 1, 2, false, adj_tris,
                                    moab::Interface::UNION);
    adj_tris = intersect(adj_tris, fineTris);
    if (adj_tris.size() != 2) {
      continue;
    }
 
    FTensor::Tensor1<double, 3> n0;
    FTensor::Tensor1<double, 3> n1;
    CHKERR tri_normal(adj_tris[0], n0);
    CHKERR tri_normal(adj_tris[1], n1);
    const double dot = n0(i) * n1(i);
    if (std::abs(dot) < fine_planar_cos) {
      fineFeatureTris.insert(adj_tris[0]);
      fineFeatureTris.insert(adj_tris[1]);
    }
  }
 
  filterOuterTrisByEdges(fineMoab, fineFeatureTris);
  if (fineFeatureTris.empty()) {
    fineFeatureTris = fineTris;
  }
  if (debugSurfaces) {
    CHKERR CutMeshInterface::SaveData(fineMoab)("fine_tri_features.vtk",
                                                fineFeatureTris);
  }
 
  // KD tree on fine surface for closest-point queries.
  fineRoot = 0;
  CHKERR fineTree.build_tree(fineTris, &fineRoot);
  fineFeatureRoot = 0;
  CHKERR fineFeatureTree.build_tree(fineFeatureTris, &fineFeatureRoot);
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode Projection10NodeCoordsFromFineSurfaceOnField::operator()() {
  MoFEMFunctionBeginHot;
 
  if (entPtr == NULL) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
  }
  if (entPtr->getName() != fieldName)
    MoFEMFunctionReturnHot(0);
  if (entPtr->getEntType() != MBEDGE)
    MoFEMFunctionReturnHot(0);
 
  const EntityHandle edge = entPtr->getEnt();
  if (coarseSkinEdges.find(edge) == coarseSkinEdges.end()) {
    MoFEMFunctionReturnHot(0);
  }
 
  Range adj_tris;
  CHKERR mField.get_moab().get_adjacencies(&edge, 1, 2, false, adj_tris,
                                          moab::Interface::UNION);
  adj_tris = intersect(adj_tris, coarseSkinTris);
  if (adj_tris.size() != 2) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "skin edge does not have exactly two adjacent triangles");
  }
 
  int num_nodes = 0;
  const EntityHandle *conn = nullptr;
  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
  if (num_nodes < 2) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "edge connectivity has fewer than 2 nodes");
  }
 
  const auto base = entPtr->getApproxBase();
  double edge_shape_function_val = 0.0;
  switch (base) {
  case AINSWORTH_LEGENDRE_BASE:
    edge_shape_function_val = 0.25;
    break;
  case AINSWORTH_LOBATTO_BASE:
    edge_shape_function_val = 0.25 * LOBATTO_PHI0(0);
    break;
  case DEMKOWICZ_JACOBI_BASE: {
    double L[3];
    CHKERR Legendre_polynomials(2, 0, NULL, L, NULL, 1);
    edge_shape_function_val = 0.125 * LOBATTO_PHI0(0);
  } break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
  }
 
  auto ent_data = entPtr->getEntFieldData();
  auto t_ent_data = getFTensor1FromPtr<3>(&ent_data(0));
  if (ent_data.size() != 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "edge data size mismatch");
  }
  coords.resize(num_nodes * 3);
  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
  FTensor::Tensor1<double, 3> edge_vec;
  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_n0(
      &coords[0], &coords[1], &coords[2]);
  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_n1(
      &coords[3], &coords[4], &coords[5]);
  edge_vec(i) = t_n1(i) - t_n0(i);
  const double h = edge_vec.l2();
  if (h <= 0.0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "zero-length edge");
  }
  aveMidCoord.resize(3);
  midNodeCoord.resize(3);
  auto t_ave_mid_coord = getFTensor1FromPtr<3>(&aveMidCoord(0));
  auto t_mid_node_coord = getFTensor1FromPtr<3>(&midNodeCoord(0));
  t_ave_mid_coord(i) = 0.5 * (t_n0(i) + t_n1(i));
  t_mid_node_coord(i) =
      t_ave_mid_coord(i) + edge_shape_function_val * t_ent_data(i);
 
  auto tri_normal = [&](const EntityHandle tri,
                        FTensor::Tensor1<double, 3> &t_normal) {
    MoFEMFunctionBeginHot;
    CHKERR mField.getInterface<Tools>()->getTriNormal(tri, &t_normal(0));
    t_normal.normalize();
    MoFEMFunctionReturnHot(0);
  };
 
  FTensor::Tensor1<double, 3> n0;
  FTensor::Tensor1<double, 3> n1;
  CHKERR tri_normal(adj_tris[0], n0);
  CHKERR tri_normal(adj_tris[1], n1);
  const double dot = n0(i) * n1(i);
 
  if (std::abs(dot) >= planarAngleCos) {
    MoFEMFunctionReturnHot(0);
  }
 
  double closest_point[3] = {0.0, 0.0, 0.0};
  EntityHandle closest_tri = 0;
  const double query[3] = {midNodeCoord[0], midNodeCoord[1], midNodeCoord[2]};
  CHKERR fineFeatureTree.closest_triangle(fineFeatureRoot, query, closest_point,
                                          closest_tri);
  if (!closest_tri) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no closest triangle found on fine mesh");
  }
 
  auto get_disp = [&](double _closest_point[3]) {
    FTensor::Tensor1<double, 3> delta;
    for (int ii = 0; ii < 3; ii++)
      delta(ii) = _closest_point[ii] - midNodeCoord[ii];
    return delta.l2();
  };
 
  double disp = get_disp(closest_point);
 
  if (disp > maxDispFactor * h) {
    // std::cout << "Fallback to full fine tree for edge " << edge << std::endl;
    EntityHandle feature_tri = 0;
    CHKERR fineTree.closest_triangle(fineRoot, query, closest_point,
                                     feature_tri);
    if (feature_tri) {
      disp = get_disp(closest_point);
    }
  }
 
  if (disp > maxDispFactor * h) {
    MOFEM_LOG("SELF", Sev::warning) <<
        "The displacement for edge " << edge
        << " is too large: " << disp << " (max allowed (-max_midnode_disp_factor * h) "
        << maxDispFactor * h << "). Using original mid-edge position.";
    
  } else {
    midNodeCoord[0] = closest_point[0];
    midNodeCoord[1] = closest_point[1];
    midNodeCoord[2] = closest_point[2];
  }
 
  diffNodeCoord.resize(3);
  auto t_diff_node_coord = getFTensor1FromPtr<3>(&diffNodeCoord(0));
  t_diff_node_coord(i) = t_mid_node_coord(i) - t_ave_mid_coord(i);
  dOf.resize(3);
  auto t_dof = getFTensor1FromPtr<3>(&dOf(0));
  t_dof(i) = t_diff_node_coord(i) / edge_shape_function_val;
  if (entPtr->getNbOfCoeffs() > 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only fields which are rank 3 or lower");
  }
  ent_data[0] = dOf[0];
  ent_data[1] = dOf[1];
  ent_data[2] = dOf[2];
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode Projection10NodeCoordsFromFineSurfaceOnField::postProcess() {
  MoFEMFunctionBeginHot;
  MoFEMFunctionReturnHot(0);
}
 
ProjectionFieldOn10NodeTet::ProjectionFieldOn10NodeTet(Interface &m_field,
                                                       std::string _fieldName,
                                                       bool set_nodes,
                                                       bool on_coords,
                                                       std::string on_tag)
    : Projection10NodeCoordsOnField(m_field, _fieldName), setNodes(set_nodes),
      onCoords(on_coords), onTag(on_tag), maxApproximationOrder(20) {}
 
Tag th;
Field_multiIndex::index<FieldName_mi_tag>::type::iterator field_it;
VectorDouble L;
VectorDouble K;
 
MoFEMErrorCode ProjectionFieldOn10NodeTet::preProcess() {
  MoFEMFunctionBeginHot;
  if (!onCoords) {
    if (onTag == "NoNE") {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "tag name not specified");
    }
    field_it = fieldsPtr->get<FieldName_mi_tag>().find(fieldName);
    if (field_it == fieldsPtr->get<FieldName_mi_tag>().end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "field not found %s",
               fieldName.c_str());
    }
    int field_rank = (*field_it)->getNbOfCoeffs();
    VectorDouble def_VAL = ublas::zero_vector<double>(field_rank);
    CHKERR mField.get_moab().tag_get_handle(
        onTag.c_str(), field_rank, MB_TYPE_DOUBLE, th,
        MB_TAG_CREAT | MB_TAG_SPARSE, &*def_VAL.data().begin());
    MOAB_THROW(rval);
  }
 
  L.resize(maxApproximationOrder + 1);
  CHKERR Legendre_polynomials(maxApproximationOrder, 0., NULL,
                              &*L.data().begin(), NULL, 3);
  K.resize(10);
  CHKERR LobattoKernel_polynomials(9, 0., NULL, &*K.data().begin(), NULL, 3);
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode ProjectionFieldOn10NodeTet::operator()() {
  MoFEMFunctionBeginHot;
  if (dofPtr->getName() != fieldName)
    MoFEMFunctionReturnHot(0);
  if (setNodes) {
    if (dofPtr->getEntType() == MBVERTEX) {
      EntityHandle node = dofPtr->getEnt();
      if (onCoords) {
        coords.resize(3);
        CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
        coords[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
        CHKERR mField.get_moab().set_coords(&node, 1, &*coords.data().begin());
      } else {
        int field_rank = (*field_it)->getNbOfCoeffs();
        if (field_rank != dofPtr->getNbOfCoeffs()) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "data inconsistency");
        }
        double *tag_value;
        int tag_size;
        CHKERR mField.get_moab().tag_get_by_ptr(
            th, &node, 1, (const void **)&tag_value, &tag_size);
        if (tag_size != dofPtr->getNbOfCoeffs()) {
          SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
        }
        tag_value[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
      }
    }
    MoFEMFunctionReturnHot(0);
  }
  if (dofPtr->getEntType() != MBEDGE) {
    MoFEMFunctionReturnHot(0);
  }
  EntityHandle edge = dofPtr->getEnt();
  if (type_from_handle(edge) != MBEDGE) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only elements which are type of MBEDGE");
  }
 
  int num_nodes;
  const EntityHandle *conn;
  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
  if ((num_nodes != 2) && (num_nodes != 3)) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "this method works only 4 node and 10 node tets");
  }
  if (num_nodes == 2) {
    MoFEMFunctionReturnHot(0);
  }
 
  if (dofPtr->getDofOrder() >= maxApproximationOrder) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "too big approximation order, increase constant "
            "max_ApproximationOrder");
  }
  double approx_val = 0;
  FieldApproximationBase base = dofPtr->getApproxBase();
  switch (base) {
  case AINSWORTH_LEGENDRE_BASE:
    approx_val = 0.25 * L[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
    break;
  case AINSWORTH_LOBATTO_BASE:
    approx_val = 0.25 * K[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
  }
 
  if (onCoords) {
    coords.resize(num_nodes * 3);
    CHKERR mField.get_moab().get_coords(conn, num_nodes,
                                        &*coords.data().begin());
    if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
      // add only one when higher order terms present
      double ave_mid = (coords[3 * 0 + dofPtr->getDofCoeffIdx()] +
                        coords[3 * 1 + dofPtr->getDofCoeffIdx()]) *
                       0.5;
      coords[2 * 3 + dofPtr->getDofCoeffIdx()] = ave_mid;
    }
    coords[2 * 3 + dofPtr->getDofCoeffIdx()] += approx_val;
    CHKERR mField.get_moab().set_coords(&conn[2], 1, &coords[3 * 2]);
  } else {
    int tag_size;
    double *tag_value[num_nodes];
    CHKERR mField.get_moab().tag_get_by_ptr(
        th, conn, num_nodes, (const void **)tag_value, &tag_size);
    if (tag_size != dofPtr->getNbOfCoeffs()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
    }
    if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
      // add only one when higher order terms present
      double ave_mid = (tag_value[0][dofPtr->getDofCoeffIdx()] +
                        tag_value[1][dofPtr->getDofCoeffIdx()]) *
                       0.5;
      tag_value[2][dofPtr->getDofCoeffIdx()] = ave_mid;
    }
    tag_value[2][dofPtr->getDofCoeffIdx()] += approx_val;
  }
  MoFEMFunctionReturnHot(0);
}
}; // namespace MoFEM