mesh_tag_transfer.cpp

Go to the documentation of this file.

/** \file mesh_tag_transfer.cpp
 
  \brief extract mesh tag
    tool to extract a mesh tag and write to new mesh file. 
    Can be used to extract a component of an existing vector tag.
 
    NOTE: Only supports extracting one component from tags defined on tetrahedral entities 
    and creates new tag on same entities. Can be extended to support other entity types as needed.
*/
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    char mesh_file_name[255] = "mesh.h5m";
    char mesh_out_file[255] = "out.h5m";
    char tag_name[255] = "new_tag";
    char extract_tag_name[255] = "old_tag_name";
    int extract_component = 0;
 
    PetscBool flg_extract_tag_name = PETSC_FALSE;
    int atom_test = 0;
 
    PetscBool flg_file = PETSC_FALSE;
    PetscBool flg_analytical = PETSC_FALSE;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
    CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_file);
    CHKERR PetscOptionsString("-create_tag_name", "name of tag to create", "",
                              "", tag_name, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-extract_tag_name",
                              "name of tag to extract from", "", "",
                              extract_tag_name, 255, &flg_extract_tag_name);
    CHKERR PetscOptionsInt("-component_to_extract",
                           "component of vector tag to extract (0-based index)",
                           "", extract_component, &extract_component,
                           PETSC_NULLPTR);
    CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
                              "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-atom_test", "flag to create tag for test", "",
                           atom_test, &atom_test, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-use_analytical_tag",
                            "bool to use coded analytical function", "",
                            flg_analytical, &flg_analytical, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);
 
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    // Create MoFEM database
    MoFEM::Core core(moab);
 
    Range ents;
    std::vector<double> new_tag_values;
    // default entity type is tets, can be updated as needed
    EntityType entity_type = MBTET;
 
    if (atom_test || flg_analytical)
      entity_type = MBVERTEX;
 
    if (flg_extract_tag_name) {
      // Tag extraction only supported on tets for now, can be extended to other entity types as needed
      CHKERR moab.get_entities_by_type(0, entity_type, ents);
      // find existing tag to extract value from
      Tag extract_tag;
      auto rval = moab.tag_get_handle(extract_tag_name, extract_tag);
      if (rval != MB_SUCCESS) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Failed to get handle for tag '%s'", extract_tag_name);
      }
 
      int extract_tag_len;
      CHKERR moab.tag_get_length(extract_tag, extract_tag_len);
      MOFEM_LOG("WORLD", Sev::inform)
          << "Extracting component " << extract_component << " from tag '"
          << extract_tag_name << "' with length " << extract_tag_len;
 
      std::vector<double> extract_values(extract_tag_len * ents.size(), 0.0);
      CHKERR moab.tag_get_data(extract_tag, ents, &extract_values[0]);
 
      MOFEM_LOG("WORLD", Sev::inform)
          << "Extracted " << extract_values.size() << " values from tag '"
          << extract_tag_name << "'";
 
      // extract specified component into new tag values
      for (size_t i = 0; i < ents.size(); ++i) {
        int idx = extract_component; // Assuming 1 component
        if (idx < extract_tag_len) {
          new_tag_values.push_back(extract_values[i * extract_tag_len + idx]);
        } else {
          SETERRQ(
              PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Component index to extract is out of bounds for the tag length");
        }
      }
      MOFEM_LOG("WORLD", Sev::inform) << "Extracted component values for "
                                      << new_tag_values.size() << " entities";
    } else {
      MOFEM_LOG("WORLD", Sev::warning) << "No tag specified for extraction.";
      std::vector<double> test_data;
      if (atom_test) {
        // Atom test assigns values based on vector index
        CHKERR moab.get_entities_by_type(0, entity_type, ents);
 
        std::vector<double> test_data;
        if (atom_test == 1) {
          for (size_t i = 0; i < ents.size(); ++i) {
            // assign index 0-5 for each
            for (int idx = 0; idx < 5; ++idx) {
              test_data.push_back(idx);
            }
          }
        } else {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "atom test '%i' does not exist for setting a tag", atom_test);
        }
 
        for (const auto &val : test_data) {
          new_tag_values.push_back(val);
        }
      } else if (flg_analytical) {
        // create new tag for youngs modulus and set value based on analytical function
        auto analytical_elastic = [](double, double,
                                     const std::vector<double> &x,
                                     const std::vector<double> &y,
                                     const std::vector<double> &z, int) {
          std::vector<double> youngs_modulus;
          youngs_modulus.reserve(x.size());
 
          for (double xi : x) {
            double E0 = 1000.0;
            double local = E0 * std::exp(xi);
            youngs_modulus.push_back(local);
          }
 
          return youngs_modulus;
        };
 
        std::vector<double> coords(3 * ents.size(), 0.0);
        CHKERR moab.get_coords(ents, &coords[0]);
 
        for (size_t i = 0; i < ents.size(); ++i) {
          double x = coords[3 * i + 0];
          double y = coords[3 * i + 1];
          double z = coords[3 * i + 2];
 
          double val = analytical_elastic(0.0, 0.0, {x}, {y}, {z}, 0)[0];
          test_data.push_back(val);
        }
 
        for (const auto &val : test_data) {
          new_tag_values.push_back(val);
        }
      } else {
        MOFEM_LOG("WORLD", Sev::warning) << "No tag has been created.";
      }
    }
 
    // print the range of new tag values for verification
    auto minmax =
        std::minmax_element(new_tag_values.begin(), new_tag_values.end());
    MOFEM_LOG("WORLD", Sev::inform)
        << "New tag '" << tag_name << "' value range: [" << *minmax.first
        << ", " << *minmax.second << "]";
 
    // Create new tag and set value
    std::vector<double> def_val;
 
    if (atom_test == 1) {
      def_val = {0.0, 0.0, 0.0, 0.0, 0.0};
    } else {
      def_val = {0.0};
    }
 
    Tag new_tag;
    CHKERR moab.tag_get_handle(tag_name, def_val.size(), MB_TYPE_DOUBLE, new_tag,
                               MB_TAG_CREAT | MB_TAG_DENSE, &def_val);
 
    CHKERR moab.tag_set_data(new_tag, ents, &new_tag_values[0]);
 
    // write output mesh with new tag
    CHKERR moab.write_file(mesh_out_file, "MOAB", "PARALLEL=WRITE_PART");
    CHKERR moab.write_file("out.vtk", "vtk");
 
    // check results
    if (atom_test == 1) {
      std::vector<double> check_values(new_tag_values.size(), 0.0);
      CHKERR moab.tag_get_data(new_tag, ents, &check_values[0]);
      MOFEM_LOG("WORLD", Sev::inform)
          << "check_values size: " << check_values.size()
          << ", new_tag_values size: " << new_tag_values.size();
 
      for (size_t i = 0; i < new_tag_values.size(); ++i) {
        if (std::abs(check_values[i] - new_tag_values[i]) > 1e-6) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Atom test failed: Value mismatch at index %zu: expected %f, "
                  "got %f",
                  i, new_tag_values[i], check_values[i]);
        }
      }
    } else if (atom_test == 2) {
      // For index test, check that extracted values match the expected component index
      // instead of new_tag_values use expected values based on component index
      std::vector<double> expected_values(new_tag_values.size(), 0.0);
      CHKERR moab.tag_get_data(new_tag, ents, &expected_values[0]);
      MOFEM_LOG("WORLD", Sev::inform)
          << "check_values size: " << expected_values.size()
          << ", new_tag_values size: " << new_tag_values.size();
 
      for (size_t i = 0; i < expected_values.size(); ++i) {
        if (std::abs(expected_values[i] - extract_component) > 1e-6) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Tag extraction test failed: Value mismatch at index %zu: "
                  "expected %i, got %f",
                  i, extract_component, expected_values[i]);
        }
      }
    } else if (atom_test > 0) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "atom test %i does not exist", atom_test);
    }
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}