constFieldAccessor.hpp

#ifndef ABLATELIBRARY_CONSTFIELDACCESSOR_HPP
#define ABLATELIBRARY_CONSTFIELDACCESSOR_HPP
 
#include <petsc.h>
#include "field.hpp"
 
namespace ablate::domain {
template <class DataType, bool isLocal = true>
class ConstFieldAccessor {
   private:
    Vec dataVector;
 
    const Field& dataField;
 
    const PetscScalar* dataArray{};
 
    DM dataDM;
 
   public:
    ConstFieldAccessor(Vec dataVectorIn, const Field& dataFieldIn, DM dm = nullptr) : dataVector(dataVectorIn), dataField(dataFieldIn), dataDM(dm) {
        // Get read/write access to the vector
        VecGetArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
 
        // Get the dm from the vector if not specified
        if (!dataDM) {
            VecGetDM(dataVector, &dataDM) >> ablate::utilities::PetscUtilities::checkError;
        }
    }
 
    ~ConstFieldAccessor() {
        // Put back the vector
        VecRestoreArrayRead(dataVector, &dataArray) >> ablate::utilities::PetscUtilities::checkError;
    };
 
    [[nodiscard]] inline const Field& GetField() const { return dataField; }
 
    template <class IndexType>
    inline const DataType* operator[](IndexType point) const {
        const DataType* field;
        if constexpr (isLocal) {
            DMPlexPointLocalFieldRead(dataDM, point, dataField.id, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
        } else {
            DMPlexPointGlobalFieldRead(dataDM, point, dataField.id, dataArray, &field) >> ablate::utilities::PetscUtilities::checkError;
        }
        return field;
    }
 
    template <class IndexType>
    inline PetscErrorCode operator()(IndexType point, const DataType** field) const {
        PetscFunctionBeginHot;
        if constexpr (isLocal) {
            PetscCall(DMPlexPointLocalFieldRead(dataDM, point, dataField.id, dataArray, field));
        } else {
            PetscCall(DMPlexPointGlobalFieldRead(dataDM, point, dataField.id, dataArray, field));
        }
        PetscFunctionReturn(PETSC_SUCCESS);
    }
 
    ConstFieldAccessor(const ConstFieldAccessor&) = delete;
};
}  // namespace ablate::domain
#endif  // ABLATELIBRARY_CONSTFIELDACCESSOR_HPP