MatOps.cpp

Go to the documentation of this file.

/**
 * @file MatOps.cpp
 * @author your name (you@domain.com)
 * @brief 
 * @date 2026-03-26
 * 
 * @copyright Copyright (c) 2026
 * 
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
#include "MatOps.hpp"
 
namespace MatOps {
 
static void sortStateTags(StateTags &state_tags) {
  std::sort(
      state_tags.begin(), state_tags.end(),
      [](const auto &lhs, const auto &rhs) { return lhs.name < rhs.name; });
}
 
int MatOpsTagsRegistry::atTagName(std::string name) {
  return tagNameVsTag.at(name);
}
 
int MatOpsTagsRegistry::setTagName(std::string name, int tag) {
#ifndef NDEBUG
  MOFEM_LOG_CHANNEL("WORLD");
#endif
 
  auto it = tagNameVsTag.find(name);
 
  if (tag >= 0) {
    if (it != tagNameVsTag.end()) {
#ifndef NDEBUG
      MOFEM_LOG("WORLD", Sev::warning)
          << "Tag name: " << name
          << " is already registered with tag: " << it->second
          << ". Overwriting with new tag: " << tag;
#endif
      it->second = -1; // set it to invalid value to avoid confusion
    }
#ifndef NDEBUG
    for (auto &p : tagNameVsTag) {
      if (p.second == tag) {
        CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                          "Tag: " + std::to_string(tag) +
                              " is already registered with name: " + p.first);
      }
    }
#endif
    // register new tag
    tagNameVsTag[name] = tag;
    return tag;
  }
 
  int max_tag = 0;
  for (auto &p : tagNameVsTag) {
    if (p.second > max_tag) {
      max_tag = p.second;
    }
  }
 
  if (it == tagNameVsTag.end()) {
    tagNameVsTag[name] = max_tag + 1;
  }
 
#ifndef NDEBUG
  MOFEM_LOG("WORLD", Sev::inform)
      << "Tag name: " << name
      << " is registered with tag: " << tagNameVsTag[name];
#endif
 
  return tagNameVsTag[name];
}
 
int MatOpsTagsRegistry::getTagByName(std::string name) {
  auto it = tagNameVsTag.find(name);
  if (it == tagNameVsTag.end()) {
    int new_tag = tagNameVsTag.size() + 1;
    tagNameVsTag[name] = new_tag;
    return new_tag;
  }
  return it->second;
}
 
std::string MatOpsTagsRegistry::getTagName(int tag) {
  for (const auto &p : tagNameVsTag) {
    if (p.second == tag) {
      return p.first;
    }
  }
  return std::string();
}
 
struct MatOpsDataImpl : public MatOpsData {
  using MatrixPtr = MatOpsData::MatrixPtr;
  using StateMatrixPtr = MatOpsData::StateMatrixPtr;
 
  boost::shared_ptr<UserData> getUserDataPtr() const override {
    return userDataPtr;
  }
 
  void setUserDataPtr(boost::shared_ptr<UserData> user_data_ptr) override {
    userDataPtr = user_data_ptr;
  }
 
  struct Data {
    Data() = default;
    Data(const std::string &data_name,
         const MatrixDouble &data_value = MatrixDouble(), int data_shift = 0)
        : name(data_name), data(data_value), shift(data_shift) {}
    Data(const std::string &data_name, int num_components, Tag data_tag = 0)
        : name(data_name), data(num_components, 1),
          numComponents(num_components), tag(data_tag) {}
 
    std::string name;
    mutable MatrixDouble data;
    mutable int shift = 0;
    int numComponents = 0;
    Tag tag = 0;
    mutable EntityHandle cachedEntity = 0;
    mutable double *cachedPtr = nullptr;
    mutable int cachedSize = 0;
 
    void resetStateCache() const {
      cachedEntity = 0;
      cachedPtr = nullptr;
      cachedSize = 0;
    }
  };
 
  struct DataSequence_mi_tag;
  struct DataName_mi_tag;
 
  using DataContainer = boost::multi_index::multi_index_container<
 
      Data,
 
      boost::multi_index::indexed_by<
 
          boost::multi_index::sequenced<
              boost::multi_index::tag<DataSequence_mi_tag>>,
 
          boost::multi_index::hashed_unique<
              boost::multi_index::tag<DataName_mi_tag>,
              boost::multi_index::member<Data, std::string, &Data::name>>>>;
 
  auto &getDataCommonByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataCommon);
  }
 
  auto &getDataActiveByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataActive);
  }
 
  auto &getDataDependentByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataDependent);
  }
 
  auto &getDataDependentDerivativesByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataDependentDerivatives);
  }
 
  auto &getStateDataByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataState);
  }
 
  auto &getStateDataBySequence() {
    return boost::multi_index::get<DataSequence_mi_tag>(dataState);
  }
 
  std::pair<bool, MatrixPtr>
  insertCommonData(const std::string &name,
                   const MatrixDouble &data = MatrixDouble(),
                   int shift = 0) override {
    auto insert_ret = getDataCommonByName().insert(Data{name, data, shift});
    return {insert_ret.second,
            boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                            &insert_ret.first->data)};
  }
 
  std::pair<bool, MatrixPtr>
  insertActiveData(const std::string &name,
                   const MatrixDouble &data = MatrixDouble(),
                   int shift = 0) override {
    auto insert_ret = getDataActiveByName().insert(Data{name, data, shift});
    return {insert_ret.second,
            boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                            &insert_ret.first->data)};
  }
 
  std::pair<bool, MatrixPtr>
  insertDependentData(const std::string &name,
                      const MatrixDouble &data = MatrixDouble(),
                      int shift = 0) override {
    auto insert_ret = getDataDependentByName().insert(Data{name, data, shift});
    return {insert_ret.second,
            boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                            &insert_ret.first->data)};
  }
 
  std::pair<bool, MatrixPtr>
  insertDependentDerivativesData(const std::string &name,
                                 const MatrixDouble &data = MatrixDouble(),
                                 int shift = 0) override {
    auto insert_ret =
        getDataDependentDerivativesByName().insert(Data{name, data, shift});
    return {insert_ret.second,
            boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                            &insert_ret.first->data)};
  }
 
  MatrixPtr getCommonDataPtr(const std::string &name) override {
    auto &data_by_name = getDataCommonByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Common data with name " + name + " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  MatrixPtr getActiveDataPtr(const std::string &name) override {
    auto &data_by_name = getDataActiveByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Active data with name " + name + " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  MatrixPtr getDependentDataPtr(const std::string &name) override {
    auto &data_by_name = getDataDependentByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Dependent data with name " + name + " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  MatrixPtr getDependentDerivativesDataPtr(const std::string &name) override {
    auto &data_by_name = getDataDependentDerivativesByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Dependent derivatives data with name " + name +
                            " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  StateMatrixPtr getStateDataPtr(const std::string &name, EntityHandle ent,
                                 int gg) override {
    auto &state_by_name = getStateDataByName();
    auto it = state_by_name.find(name);
    if (it == state_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "State data with name " + name + " not found");
      return StateMatrixPtr();
    }
    if (!mFieldPtr) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "MoFEM interface for state data is not set");
      return StateMatrixPtr();
    }
    if (!ent) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "State entity is not set");
      return StateMatrixPtr();
    }
    if (gg < 0) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "State gauss point out of range");
      return StateMatrixPtr();
    }
    if (!it->tag) {
      CHK_THROW_MESSAGE(
          MOFEM_OPERATION_UNSUCCESSFUL,
          "State data " + name +
              " is not initialised; call setupStateData() after binding");
      return StateMatrixPtr();
    }
    const auto required_size = it->numComponents * (gg + 1);
    if (gg == 0 || ent != it->cachedEntity || !it->cachedPtr ||
        required_size > it->cachedSize) {
      it->resetStateCache();
      double *tag_data = nullptr;
      int tag_size = 0;
      auto &moab = mFieldPtr->get_moab();
      auto rval = moab.tag_get_by_ptr(it->tag, &ent, 1,
                                      (const void **)&tag_data, &tag_size);
      if (tag_data && tag_size > 0 && tag_size % it->numComponents != 0) {
        CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                          "Inconsistent state tag size for " + name);
        return StateMatrixPtr();
      }
      if (rval != MB_SUCCESS || !tag_data || tag_size < required_size) {
        std::vector<double> storage(required_size, 0);
        if (rval == MB_SUCCESS && tag_data && tag_size > 0) {
          std::copy_n(tag_data, std::min(tag_size, required_size),
                      storage.data());
        }
        void const *storage_ptr[] = {storage.data()};
        int storage_size = required_size;
        CHK_THROW_MESSAGE(
            moab.tag_set_by_ptr(it->tag, &ent, 1, storage_ptr, &storage_size),
            "Failed to set state data " + name);
        CHK_THROW_MESSAGE(moab.tag_get_by_ptr(it->tag, &ent, 1,
                                              (const void **)&tag_data,
                                              &tag_size),
                          "Failed to get state data " + name);
      }
      if (!tag_data || tag_size < required_size ||
          (tag_size > 0 && tag_size % it->numComponents != 0)) {
        CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                          "Inconsistent state tag size for " + name);
        return StateMatrixPtr();
      }
 
      it->cachedEntity = ent;
      it->cachedPtr = tag_data;
      it->cachedSize = tag_size;
    }
 
    auto *ptr = it->cachedPtr + gg * it->numComponents;
    return boost::make_shared<MatrixAdaptor>(
        it->numComponents, 1,
        ublas::shallow_array_adaptor<double>(it->numComponents, ptr));
  }
 
  StateTags getStateTags() const override {
    StateTags state_tags;
    const auto &state_by_sequence =
        boost::multi_index::get<DataSequence_mi_tag>(dataState);
 
    for (const auto &state_data : state_by_sequence) {
      if (!state_data.tag || state_data.name.empty() ||
          state_data.numComponents <= 0)
        continue;
      state_tags.push_back(
          {state_data.name, state_data.tag, state_data.numComponents});
    }
 
    sortStateTags(state_tags);
    return state_tags;
  }
 
  MoFEMErrorCode bindStateTag(MoFEM::Interface &m_field,
                              const std::string &name,
                              int num_components) override {
    MoFEMFunctionBegin;
    for (const auto &state_data : getStateDataBySequence()) {
      state_data.resetStateCache();
    }
    mFieldPtr = &m_field;
 
    if (name.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "State data name can not be empty");
    }
 
    if (num_components <= 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "State data %s has non-positive number of components %d",
              name.c_str(), num_components);
    }
 
    auto &state_by_name = getStateDataByName();
    auto it = state_by_name.find(name);
    if (it != state_by_name.end()) {
      if (it->numComponents != num_components) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "State data %s already exists with incompatible size",
                name.c_str());
      }
    } else {
      state_by_name.insert(Data{name, num_components});
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setupStateData() override {
    MoFEMFunctionBegin;
    if (!mFieldPtr) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "MoFEM interface for state data is not set");
    }
 
    for (auto &state_data : getStateDataBySequence()) {
      state_data.resetStateCache();
      Tag tag = 0;
      const auto state_tag_name = "_ADOLC_STATE_" +
                                  std::to_string(state_data.numComponents) +
                                  "_" + state_data.name;
      CHKERR mFieldPtr->get_moab().tag_get_handle(
          state_tag_name.c_str(), 0, MB_TYPE_DOUBLE, tag,
          MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, PETSC_NULLPTR);
      getStateDataBySequence().modify(
          getStateDataBySequence().iterator_to(state_data),
          [&](auto &data) { data.tag = tag; });
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode
  setActiveContinuousVector(std::vector<double> &active_variables) override {
    return setContinuousVector(dataActive, active_variables);
  }
 
  MoFEMErrorCode setDependentContinuousVector(
      std::vector<double> &dependent_variables) override {
    return setContinuousVector(dataDependent, dependent_variables);
  }
 
  MoFEMErrorCode setDependentDerivativesContinuousVector(
      std::vector<double> &dependent_variables_derivatives) override {
    return setContinuousVector(dataDependentDerivatives,
                               dependent_variables_derivatives);
  }
 
  MoFEMErrorCode getActiveContinuousVector(
      const std::vector<double> &active_variables) override {
    return getContinuousVector(dataActive, active_variables);
  }
 
  MoFEMErrorCode getDependentContinuousVector(
      const std::vector<double> &dependent_variables) override {
    return getContinuousVector(dataDependent, dependent_variables);
  }
 
  MoFEMErrorCode getDependentDerivativesContinuousVector(
      const std::vector<double> &dependent_variables_derivatives) override {
    return getContinuousVector(dataDependentDerivatives,
                               dependent_variables_derivatives);
  }
 
protected:
  MoFEMErrorCode setContinuousVector(DataContainer &data_container,
                                     std::vector<double> &variables) {
    MoFEMFunctionBeginHot;
 
    auto &data_by_sequence =
        boost::multi_index::get<DataSequence_mi_tag>(data_container);
 
    int nb_variables = 0;
    for (const auto &d : data_by_sequence)
      nb_variables += d.data.data().size();
 
    variables.resize(nb_variables);
 
    int shift = 0;
    for (const auto &d : data_by_sequence) {
      d.shift = shift;
      for (auto it = d.data.data().begin(); it != d.data.data().end();
           ++it, ++shift)
        variables[shift] = *it;
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
  MoFEMErrorCode getContinuousVector(DataContainer &data_container,
                                     const std::vector<double> &variables) {
    MoFEMFunctionBeginHot;
 
    auto &data_by_sequence =
        boost::multi_index::get<DataSequence_mi_tag>(data_container);
 
    int nb_variables = 0;
    for (const auto &d : data_by_sequence)
      nb_variables += d.data.data().size();
 
    if (variables.size() != static_cast<size_t>(nb_variables)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Inconsistent size of continuous vector");
    }
 
    int shift = 0;
    for (const auto &d : data_by_sequence) {
      d.shift = shift;
      for (auto it = d.data.data().begin(); it != d.data.data().end();
           ++it, ++shift)
        *it = variables[shift];
    }
 
    MoFEMFunctionReturnHot(0);
  }
 
  DataContainer dataCommon;
  DataContainer dataActive;
  DataContainer dataDependent;
  DataContainer dataDependentDerivatives;
  DataContainer dataState;
  MoFEM::Interface *mFieldPtr = nullptr;
  boost::shared_ptr<UserData> userDataPtr;
};
 
struct MatOpsDataDerivedImpl : public MatOpsDataImpl {
 
  MatOpsDataDerivedImpl(boost::weak_ptr<MatOpsData> base_data_ptr,
                        DataDictionary common_data_dictionary,
                        DataDictionary state_data_dictionary)
      : baseImplPtr(lockBaseImplPtr(base_data_ptr)),
        commonDataDictionary(common_data_dictionary),
        stateDataDictionary(state_data_dictionary) {}
 
  std::pair<bool, MatrixPtr>
  insertCommonData(const std::string &name,
                   const MatrixDouble &data = MatrixDouble(),
                   int shift = 0) override {
    return baseImplPtr->insertCommonData(name, data, shift);
  }
 
  auto &getDataActiveByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataActive);
  }
 
  auto &getDataDependentByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataDependent);
  }
 
  auto &getDataDependentDerivativesByName() {
    return boost::multi_index::get<DataName_mi_tag>(dataDependentDerivatives);
  }
 
  std::pair<bool, MatrixPtr>
  insertActiveData(const std::string &name,
                   const MatrixDouble &data = MatrixDouble(),
                   int shift = 0) override {
    auto insert_ret = getDataActiveByName().insert(Data{name, data, shift});
    if (!insert_ret.second)
      return {false, MatrixPtr()};
    return {true, boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                                  &insert_ret.first->data)};
  }
 
  std::pair<bool, MatrixPtr>
  insertDependentData(const std::string &name,
                      const MatrixDouble &data = MatrixDouble(),
                      int shift = 0) override {
    auto insert_ret = getDataDependentByName().insert(Data{name, data, shift});
    if (!insert_ret.second)
      return {false, MatrixPtr()};
    return {true, boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                                  &insert_ret.first->data)};
  }
 
  std::pair<bool, MatrixPtr>
  insertDependentDerivativesData(const std::string &name,
                                 const MatrixDouble &data = MatrixDouble(),
                                 int shift = 0) override {
    auto insert_ret =
        getDataDependentDerivativesByName().insert(Data{name, data, shift});
    if (!insert_ret.second)
      return {false, MatrixPtr()};
    return {true, boost::shared_ptr<MatrixDouble>(this->shared_from_this(),
                                                  &insert_ret.first->data)};
  }
 
  MatrixPtr getCommonDataPtr(const std::string &name) override {
    if (auto it = commonDataDictionary.find(name);
        it != commonDataDictionary.end()) {
      return baseImplPtr->getCommonDataPtr(it->second);
    }
    return baseImplPtr->getCommonDataPtr(name);
  }
 
  MatrixPtr getActiveDataPtr(const std::string &name) override {
    auto &data_by_name = getDataActiveByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Active data with name " + name + " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  MatrixPtr getDependentDataPtr(const std::string &name) override {
    auto &data_by_name = getDataDependentByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Dependent data with name " + name + " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  MatrixPtr getDependentDerivativesDataPtr(const std::string &name) override {
    auto &data_by_name = getDataDependentDerivativesByName();
    auto it = data_by_name.find(name);
    if (it == data_by_name.end()) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Dependent derivatives data with name " + name +
                            " not found");
      return MatrixPtr();
    }
    return boost::shared_ptr<MatrixDouble>(this->shared_from_this(), &it->data);
  }
 
  StateMatrixPtr getStateDataPtr(const std::string &name, EntityHandle ent,
                                 int gg) override {
    if (auto it = stateDataDictionary.find(name);
        it != stateDataDictionary.end())
      return baseImplPtr->getStateDataPtr(it->second, ent, gg);
    return baseImplPtr->getStateDataPtr(name, ent, gg);
  }
 
  StateTags getStateTags() const override {
    auto state_tags = baseImplPtr->getStateTags();
    for (auto &tag_info : state_tags) {
      for (const auto &[alias_name, base_name] : stateDataDictionary) {
        if (base_name == tag_info.name) {
          tag_info.name = alias_name;
          break;
        }
      }
    }
    sortStateTags(state_tags);
    return state_tags;
  }
 
  MoFEMErrorCode bindStateTag(MoFEM::Interface &m_field,
                              const std::string &name,
                              int num_components) override {
    if (auto it = stateDataDictionary.find(name);
        it != stateDataDictionary.end())
      return baseImplPtr->bindStateTag(m_field, it->second, num_components);
    return baseImplPtr->bindStateTag(m_field, name, num_components);
  }
 
  MoFEMErrorCode setupStateData() override {
    return baseImplPtr->setupStateData();
  }
 
  MoFEMErrorCode
  setActiveContinuousVector(std::vector<double> &active_variables) override {
    return setContinuousVector(dataActive, active_variables);
  }
 
  MoFEMErrorCode setDependentContinuousVector(
      std::vector<double> &dependent_variables) override {
    return setContinuousVector(dataDependent, dependent_variables);
  }
 
  MoFEMErrorCode setDependentDerivativesContinuousVector(
      std::vector<double> &dependent_variables_derivatives) override {
    return setContinuousVector(dataDependentDerivatives,
                               dependent_variables_derivatives);
  }
 
  MoFEMErrorCode getActiveContinuousVector(
      const std::vector<double> &active_variables) override {
    return getContinuousVector(dataActive, active_variables);
  }
 
  MoFEMErrorCode getDependentContinuousVector(
      const std::vector<double> &dependent_variables) override {
    return getContinuousVector(dataDependent, dependent_variables);
  }
 
  MoFEMErrorCode getDependentDerivativesContinuousVector(
      const std::vector<double> &dependent_variables_derivatives) override {
    return getContinuousVector(dataDependentDerivatives,
                               dependent_variables_derivatives);
  }
 
protected:
  DataContainer dataActive;
  DataContainer dataDependent;
  DataContainer dataDependentDerivatives;
 
  boost::shared_ptr<MatOpsDataImpl> baseImplPtr;
  DataDictionary commonDataDictionary;
  DataDictionary stateDataDictionary;
 
private:
  static boost::shared_ptr<MatOpsDataImpl>
  lockBaseImplPtr(boost::weak_ptr<MatOpsData> base_data_ptr) {
    auto base_data = base_data_ptr.lock();
    if (!base_data) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Base MatOps data pointer expired");
    }
 
    auto base_impl_ptr = boost::dynamic_pointer_cast<MatOpsDataImpl>(base_data);
    if (!base_impl_ptr) {
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                        "Base MatOps data pointer has incompatible type");
    }
 
    return base_impl_ptr;
  }
};
 
boost::shared_ptr<MatOpsData> createMatOpsDataPtr() {
  return boost::make_shared<MatOpsDataImpl>();
}
 
boost::shared_ptr<MatOpsData>
createMatOpsDataPtr(boost::weak_ptr<MatOpsData> base_data_ptr,
                    DataDictionary common_data_dictionary,
                    DataDictionary state_data_dictionary) {
  return boost::make_shared<MatOpsDataDerivedImpl>(
      base_data_ptr, common_data_dictionary, state_data_dictionary);
}
 
MoFEMErrorCode AdolCEvaluation::evaluateVariable(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag,
    EntityHandle entity, int gg) {
  MoFEMFunctionBegin;
  CHKERR mat_ops_data_ptr->setActiveContinuousVector(
      AdolCEvaluation::activeVariables);
  CHKERR mat_ops_data_ptr->setDependentContinuousVector(
      AdolCEvaluation::dependentVariables);
  constexpr int adolc_return_value = 0;
  int r = ::function(tag, dependentVariables.size(), activeVariables.size(),
                     activeVariables.data(), dependentVariables.data());
  if (PetscUnlikely(r < adolc_return_value)) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  CHKERR mat_ops_data_ptr->getDependentContinuousVector(
      AdolCEvaluation::dependentVariables);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode AdolCEvaluation::evaluateDerivatives(
    boost::shared_ptr<MatOpsData> mat_ops_data_ptr, int tag,
    EntityHandle entity, int gg) {
  MoFEMFunctionBegin;
  constexpr int adolc_return_value = 0;
  CHKERR mat_ops_data_ptr->setActiveContinuousVector(
      AdolCEvaluation::activeVariables);
  CHKERR mat_ops_data_ptr->setDependentContinuousVector(
      AdolCEvaluation::dependentVariables);
 
  const auto number_of_active_variables =
      AdolCEvaluation::activeVariables.size();
  const auto number_of_dependent_variables =
      AdolCEvaluation::dependentVariables.size();
  AdolCEvaluation::dependentVariablesDerivatives.resize(
      number_of_dependent_variables * number_of_active_variables);
  AdolCEvaluation::jacPtrVec.resize(number_of_dependent_variables);
  for (unsigned int n = 0; n != number_of_dependent_variables; ++n) {
    AdolCEvaluation::jacPtrVec[n] =
        &(AdolCEvaluation::dependentVariablesDerivatives
              [n * number_of_active_variables]);
  }
 
  int r = ::jacobian(tag, dependentVariables.size(), activeVariables.size(),
                     activeVariables.data(), jacPtrVec.data());
  if (PetscUnlikely(r < adolc_return_value)) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  CHKERR mat_ops_data_ptr->getDependentDerivativesContinuousVector(
      AdolCEvaluation::dependentVariablesDerivatives);
  MoFEMFunctionReturn(0);
}
 
} // namespace MatOps