subDomain.hpp

#ifndef ABLATELIBRARY_SUBDOMAIN_HPP
#define ABLATELIBRARY_SUBDOMAIN_HPP
#include <petsc.h>
#include <algorithm>
#include <map>
#include <mathFunctions/fieldFunction.hpp>
#include <memory>
#include <string>
#include "constFieldAccessor.hpp"
#include "domain.hpp"
#include "fieldAccessor.hpp"
#include "fieldDescription.hpp"
#include "io/serializable.hpp"
#include "range.hpp"
#include "utilities/petscUtilities.hpp"
 
namespace ablate::domain {
 
class SubDomain : public io::Serializable {
   private:
    inline static const std::string defaultName = "domain";
 
    Domain& domain;
 
    std::string name;
 
    DMLabel label;
 
    PetscInt labelValue;
 
    IS fieldMap;
 
    std::map<std::string, Field> fieldsByName;
 
    std::map<FieldLocation, std::vector<Field>> fieldsByType;
 
    PetscDS discreteSystem;
    PetscDS auxDiscreteSystem{};
 
    DM auxDM;
 
    Vec auxLocalVec;
 
    Vec auxGlobalVec;
 
    DM subDM;
    Vec subSolutionVec;
    DM subAuxDM;
    Vec subAuxVec;
 
    std::vector<std::shared_ptr<mathFunctions::FieldFunction>> exactSolutions;
 
    void CopyGlobalToSubVector(DM subDM, DM gDM, Vec subVec, Vec globVec, const std::vector<Field>& subFields, const std::vector<Field>& gFields = {}, bool localVector = false) const;
 
    void CopySubVectorToGlobal(DM subDM, DM gDM, Vec subVec, Vec globVec, const std::vector<Field>& subFields, const std::vector<Field>& gFields = {}, bool localVector = false) const;
 
   public:
    SubDomain(Domain& domain, PetscInt dsNumber, const std::vector<std::shared_ptr<FieldDescription>>& allAuxFields);
    ~SubDomain() override;
 
    // prevent unintended copy of the subdomain
    SubDomain(const SubDomain& temp_obj) = delete;
 
    // prevent unintended copy of the subdomain
    SubDomain& operator=(const SubDomain& temp_obj) = delete;
 
    void CreateSubDomainStructures();
 
    [[nodiscard]] inline const Field& GetField(const std::string& fieldName) const {
        if (fieldsByName.count(fieldName)) {
            return fieldsByName.at(fieldName);
        } else {
            throw std::invalid_argument("Cannot locate field " + fieldName + " in subDomain " + name);
        }
    }
 
    [[nodiscard]] inline const Field& GetField(PetscInt id, FieldLocation location = FieldLocation::SOL) const {
        auto field = std::find_if(fieldsByName.begin(), fieldsByName.end(), [location, id](auto pair) { return pair.second.location == location && pair.second.id == id; });
        if (field != fieldsByName.end()) {
            return field->second;
        } else {
            throw std::invalid_argument("Cannot locate field with id " + std::to_string(id) + " in subDomain " + name);
        }
    }
 
    [[nodiscard]] inline const std::vector<Field>& GetFields(FieldLocation type = FieldLocation::SOL) const { return fieldsByType.at(type); }
 
    [[nodiscard]] std::vector<Field> GetFields(FieldLocation type, std::string_view tag) const;
 
    [[nodiscard]] inline DM GetFieldDM(const Field& field) const noexcept {
        switch (field.location) {
            case FieldLocation::SOL:
                return GetDM();
            case FieldLocation::AUX:
                return GetAuxDM();
        }
        return nullptr;
    }
 
    [[nodiscard]] inline Vec GetVec(const Field& field) const noexcept {
        switch (field.location) {
            case FieldLocation::SOL:
                return GetSolutionVector();
            case FieldLocation::AUX:
                return GetAuxVector();
        }
        return nullptr;
    }
 
    inline Vec GetGlobalVec(const Field& field) noexcept {
        switch (field.location) {
            case FieldLocation::SOL:
                return GetSolutionVector();
            case FieldLocation::AUX:
                return GetAuxGlobalVector();
        }
        return nullptr;
    }
 
    [[nodiscard]] inline auto GetSolutionAccessor(const Field& field) const {
        if (field.location != FieldLocation::SOL) {
            throw std::invalid_argument("ablate::domain::SubDomain::GetSolutionAccessor requires a Solution Field");
        }
        return FieldAccessor<PetscScalar, false>(GetSolutionVector(), field, GetDM());
    }
 
    [[nodiscard]] inline auto GetConstSolutionAccessor(const Field& field) const {
        if (field.location != FieldLocation::SOL) {
            throw std::invalid_argument("ablate::domain::SubDomain::GetSolutionAccessor requires a Solution Field");
        }
        return ConstFieldAccessor<PetscScalar, false>(GetSolutionVector(), field, GetDM());
    }
 
    [[nodiscard]] inline auto GetAuxAccessor(const Field& field) const {
        if (field.location != FieldLocation::AUX) {
            throw std::invalid_argument("ablate::domain::SubDomain::GetSolutionAccessor requires a Aux Field");
        }
        return FieldAccessor<PetscScalar>(GetAuxVector(), field, GetAuxDM());
    }
 
    [[nodiscard]] inline auto GetConstAuxAccessor(const Field& field) const {
        if (field.location != FieldLocation::AUX) {
            throw std::invalid_argument("ablate::domain::SubDomain::GetSolutionAccessor requires a Aux Field");
        }
        return ConstFieldAccessor<PetscScalar>(GetAuxVector(), field, GetAuxDM());
    }
 
    [[nodiscard]] inline auto GetSolutionAccessor(const std::string& fieldName) const { return GetSolutionAccessor(GetField(fieldName)); }
 
    [[nodiscard]] inline auto GetConstSolutionAccessor(const std::string& fieldName) const { return GetConstSolutionAccessor(GetField(fieldName)); }
 
    [[nodiscard]] inline auto GetAuxAccessor(const std::string& fieldName) { return GetAuxAccessor(GetField(fieldName)); }
 
    [[nodiscard]] inline auto GetConstAuxAccessor(const std::string& fieldName) { return GetAuxAccessor(GetField(fieldName)); }
 
    inline bool ContainsField(const std::string& fieldName) { return fieldsByName.count(fieldName) > 0; }
 
    [[nodiscard]] bool InRegion(const domain::Region&) const;
 
    [[nodiscard]] inline bool InRegion(PetscInt point) const {
        if (!label) {
            return true;
        }
        PetscInt ptValue;
        DMLabelGetValue(label, point, &ptValue) >> utilities::PetscUtilities::checkError;
        return ptValue == labelValue;
    }
 
    inline MPI_Comm GetComm() const { return PetscObjectComm((PetscObject)domain.GetDM()); }
 
    inline DMLabel GetLabel() { return label; }
 
    inline const std::string& GetId() const override { return name; }
 
    inline PetscInt GetDimensions() const noexcept { return domain.GetDimensions(); }
 
    PetscObject GetPetscFieldObject(const Field& field);
 
    [[nodiscard]] inline DM& GetDM() const noexcept { return domain.GetDM(); }
 
    [[nodiscard]] inline DM GetAuxDM() const noexcept { return auxDM; }
 
    [[nodiscard]] inline Vec GetSolutionVector() const noexcept { return domain.GetSolutionVector(); }
 
    [[nodiscard]] inline Vec GetAuxVector() const noexcept { return auxLocalVec; }
 
    Vec GetAuxGlobalVector();
 
    DM GetSubDM();
 
    Vec GetSubSolutionVector();
 
    DM GetSubAuxDM();
 
    Vec GetSubAuxVector();
 
    inline PetscDS GetDiscreteSystem() { return discreteSystem; }
 
    inline PetscDS GetAuxDiscreteSystem() { return auxDiscreteSystem; }
 
    [[nodiscard]] inline PetscInt GetNumberFields() const { return (PetscInt)fieldsByName.size(); }
 
    void ProjectFieldFunctionsToLocalVector(const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& fieldFunctions, Vec locVec, PetscReal time = 0.0) const;
 
    PetscErrorCode ProjectFieldFunctionsToSubDM(const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& initialization, Vec globVec, PetscReal time = 0.0);
 
    void SetsExactSolutions(const std::vector<std::shared_ptr<mathFunctions::FieldFunction>>& exactSolutions);
 
    PetscErrorCode GetFieldGlobalVector(const Field&, IS* vecIs, Vec* vec, DM* subdm);
 
    PetscErrorCode RestoreFieldGlobalVector(const Field&, IS* vecIs, Vec* vec, DM* subdm);
 
    PetscErrorCode GetFieldLocalVector(const Field&, PetscReal time, IS* vecIs, Vec* vec, DM* subdm);
 
    PetscErrorCode RestoreFieldLocalVector(const Field&, IS* vecIs, Vec* vec, DM* subdm);
 
    PetscErrorCode Save(PetscViewer viewer, PetscInt sequenceNumber, PetscReal time) override;
 
    PetscErrorCode Restore(PetscViewer viewer, PetscInt sequenceNumber, PetscReal time) override;
 
    void CreateEmptySubDM(DM* inDM, std::shared_ptr<domain::Region> region = {});
 
    void GetRange(const std::shared_ptr<ablate::domain::Region>& region, PetscInt depth, ablate::domain::Range& range) const { ablate::domain::GetRange(this->GetDM(), region, depth, range); }
 
    void GetCellRange(const std::shared_ptr<ablate::domain::Region>& region, ablate::domain::Range& cellRange) const { ablate::domain::GetCellRange(this->GetDM(), region, cellRange); }
 
    void GetFaceRange(const std::shared_ptr<ablate::domain::Region>& region, ablate::domain::Range& faceRange) const { ablate::domain::GetFaceRange(this->GetDM(), region, faceRange); }
 
    void RestoreRange(ablate::domain::Range& range) const { ablate::domain::RestoreRange(range); }
};
 
}  // namespace ablate::domain
#endif  // ABLATELIBRARY_SUBDOMAIN_HPP