test_cache_on_entities.cpp

Go to the documentation of this file.

/**
 * \file test_cache_on_entities.cpp
 * \example test_cache_on_entities.cpp
 *
 * Testing entities cache. Entities cache is data abstraction enabling user to
 * pass data between operators, finite elements, and problems, in convenient
 * way.
 *
 * It can be used to store indices, or history variables, and any other data,
 * which I can not imagine at that moment of time, and how to make use of it
 * depends on user creativity and imagination.
 *
 */
 
 
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
using VolEle = VolumeElementForcesAndSourcesCore;
using VolOp = VolEle::UserDataOperator;
using FaceEle = FaceElementForcesAndSourcesCore;
using FaceOp = FaceEle::UserDataOperator;
 
struct MyStorage : EntityStorage {
  MyStorage(VectorInt &data) : globalIndices(data) {}
  VectorInt globalIndices;
};
 
/**
 * @brief Operator set cache stored data, in this example, global indices, but
 * it can be any structure
 *
 */
struct OpVolumeSet : public VolOp {
  OpVolumeSet(const std::string &field_name) : VolOp(field_name, OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
 
    MoFEMFunctionBegin;
    MOFEM_LOG_CHANNEL("SYNC");
    MOFEM_LOG_TAG("SYNC", "OpVolumeSet");
 
    // Clear data when start process element
    if (type == MBVERTEX)
      entsIndices.clear();
 
    // Get field entity pointer
    auto field_ents = data.getFieldEntities();
    if (auto e_ptr = field_ents[0]) {
      // Add indices to global storage
      entsIndices.push_back(boost::make_shared<MyStorage>(data.getIndices()));
      // Store pointer to data on entity
      e_ptr->getWeakStoragePtr() = entsIndices.back();
 
      MOFEM_LOG("SYNC", Sev::inform) << "Set " << e_ptr->getEntTypeName() << " "
                                     << side << " : " << entsIndices.size();
 
      // Check if all works
      if (auto ptr = e_ptr->getSharedStoragePtr<EntityStorage>()) {
 
        if (auto cast_ptr = boost::dynamic_pointer_cast<MyStorage>(ptr))
          MOFEM_LOG("SYNC", Sev::noisy) << "Cast works";
        else
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cast not works");
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Pointer to stored data is expected to be set");
      }
 
      if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
        MOFEM_LOG("SYNC", Sev::verbose)
            << data.getIndices() << " : " << ptr->globalIndices;
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Pointer to stored data is expected to be set");
      }
 
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Pointer to entity should be set by finite element");
    }
 
    MoFEMFunctionReturn(0);
  }
 
  using SharedVecInt = boost::shared_ptr<MyStorage>;
  static std::vector<SharedVecInt>
      entsIndices; ///< This is global static storage
};
 
std::vector<OpVolumeSet::SharedVecInt> OpVolumeSet::entsIndices;
 
/**
 * @brief Test if cached data can be accessed, and check consistency of data
 *
 */
struct OpVolumeTest : public VolOp {
  OpVolumeTest(const std::string &field_name) : VolOp(field_name, OPROW) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
 
    MoFEMFunctionBegin;
    MOFEM_LOG_CHANNEL("SYNC");
    MOFEM_LOG_TAG("SYNC", "OpVolumeTest");
 
    // Get pointer to field entities
    auto field_ents = data.getFieldEntities();
    if (auto e_ptr = field_ents[0]) {
 
      MOFEM_LOG("SYNC", Sev::inform)
          << "Test " << e_ptr->getEntTypeName() << " " << side;
 
      // Check if data are cached on entity, and if code is correct, data should
      // accessible.
      if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
 
        MOFEM_LOG("SYNC", Sev::verbose)
            << data.getIndices() << " : " << ptr->globalIndices;
 
        // Check constancy of data. Stored data are indices, and expected stored
        // that should be indices, thus difference between should be zero.
        auto diff = data.getIndices() - ptr->globalIndices;
        auto dot = inner_prod(diff, diff);
        if (dot > 0)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Test failed");
 
      } else {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Pointer to stored data is expected to be set");
      }
 
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Pointer to entity should be set by finite element");
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
struct OpVolumeAssemble : public VolOp {
  OpVolumeAssemble(const std::string &field_name, Vec v)
      : VolOp(field_name, OPROW), V(v) {}
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      std::vector<double> v(nb_dofs, 1);
      CHKERR VecSetValues<EssentialBcStorage>(V, data, &v[0], ADD_VALUES);
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  Vec V;
};
 
struct VolRule {
  int operator()(int, int, int) const { return 1; }
};
 
int main(int argc, char *argv[]) {
 
  // initialize petsc
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    // Create MoAB database
    moab::Core moab_core;
    moab::Interface &moab = moab_core;
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Simple interface
    Simple *simple_interface;
    CHKERR m_field.getInterface(simple_interface);
    {
      // get options from command line
      CHKERR simple_interface->getOptions();
      // load mesh file
      CHKERR simple_interface->loadFile();
      // add fields
      CHKERR simple_interface->addDomainField("FIELD", H1,
                                              AINSWORTH_LEGENDRE_BASE, 3);
      // set fields order
      CHKERR simple_interface->setFieldOrder("FIELD", 4);
      // setup problem
      CHKERR simple_interface->setUp();
 
      // get dm
      auto dm = simple_interface->getDM();
      // create elements
      auto domain_fe = boost::make_shared<VolEle>(m_field);
      // set integration rule
      domain_fe->getRuleHook = VolRule();
 
      auto get_skin = [&]() {
        Range ents;
        CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, ents, true);
        Skinner skin(&m_field.get_moab());
        Range skin_ents;
        CHKERR skin.find_skin(0, ents, false, skin_ents);
        // filter not owned entities, those are not on boundary
        ParallelComm *pcomm =
            ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
        Range proc_skin;
        CHKERR pcomm->filter_pstatus(skin_ents,
                                     PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                     PSTATUS_NOT, -1, &proc_skin);
        Range adj;
        for (auto d : {0, 1, 2})
          CHKERR m_field.get_moab().get_adjacencies(proc_skin, d, false, adj,
                                                    moab::Interface::UNION);
        // proc_skin.merge(adj);
        return proc_skin;
      };
 
      auto get_mark_skin_dofs = [&](Range &&skin) {
        auto problem_manager = m_field.getInterface<ProblemsManager>();
        auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
        problem_manager->markDofs(simple_interface->getProblemName(), ROW,
                                  ProblemsManager::OR, skin,
                                  *marker_ptr);
        return marker_ptr;
      };
 
      SmartPetscObj<Vec> v;
      CHKERR DMCreateGlobalVector_MoFEM(dm, v);
      auto mark_dofs = get_mark_skin_dofs(get_skin());
 
      // set operator to the volume elements
      domain_fe->getOpPtrVector().push_back(new OpVolumeSet("FIELD"));
      domain_fe->getOpPtrVector().push_back(new OpVolumeTest("FIELD"));
 
      domain_fe->getOpPtrVector().push_back(
          new OpSetBc("FIELD", true, mark_dofs));
      domain_fe->getOpPtrVector().push_back(new OpVolumeAssemble("FIELD", v));
      domain_fe->getOpPtrVector().push_back(new OpUnSetBc("FIELD"));
 
      // make integration in volume (here real calculations starts)
      CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                      domain_fe);
      OpVolumeSet::entsIndices.clear();
 
      MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
      CHKERR VecAssemblyBegin(v);
      CHKERR VecAssemblyEnd(v);
      CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
 
      CHKERR DMoFEMMeshToLocalVector(simple_interface->getDM(), v,
                                     INSERT_VALUES, SCATTER_REVERSE);
 
      // CHKERR VecView(v,PETSC_VIEWER_STDOUT_WORLD);
 
      const MoFEM::Problem *problem_ptr;
      CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
 
      double *array;
      CHKERR VecGetArray(v, &array);
      for (auto dof : *(problem_ptr->getNumeredRowDofsPtr())) {
        auto loc_idx = dof->getPetscLocalDofIdx();
 
        if (dof->getFieldData() != array[loc_idx])
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Data inconsistency");
 
        if ((*mark_dofs)[loc_idx]) {
          if (array[loc_idx] != 0)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Dof assembled, but it should be not");
        } else {
          if (array[loc_idx] == 0)
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Dof not assembled, but it should be");
        }
      }
 
      CHKERR VecRestoreArray(v, &array);
    }
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}