regionManager.h

/******************************************************************************
 *   Copyright (C) 2008-2024 by the GIMLi development team                    *
 *   Carsten Rücker carsten@resistivity.net                                   *
 *   Thomas Günther thomas@resistivity.net                                    *
 *                                                                            *
 *   Licensed under the Apache License, Version 2.0 (the "License");          *
 *   you may not use this file except in compliance with the License.         *
 *   You may obtain a copy of the License at                                  *
 *                                                                            *
 *       http://www.apache.org/licenses/LICENSE-2.0                           *
 *                                                                            *
 *   Unless required by applicable law or agreed to in writing, software      *
 *   distributed under the License is distributed on an "AS IS" BASIS,        *
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
 *   See the License for the specific language governing permissions and      *
 *   limitations under the License.                                           *
 *                                                                            *
 ******************************************************************************/
 
#ifndef _GIMLI_REGIONMANAGER__H
#define _GIMLI_REGIONMANAGER__H
 
#include "gimli.h"
#include "vector.h"
#include "trans.h"
 
#include <set>
#include <list>
 
namespace GIMLI{
 
class PosFunctor{
public:
    virtual ~PosFunctor(){}
    virtual const double operator()(const RVector3 & pos) const = 0;
};
 
class DLLEXPORT Region{
public:
    Region(SIndex marker, RegionManager * parent, bool single=false);
 
    Region(SIndex marker, const Mesh & mesh, RegionManager * parent);
 
    Region(SIndex marker, const Mesh & mesh, SIndex cellMarker, RegionManager * parent);
 
    Region(const Region & region);
 
    Region & operator = (const Region & region);
 
    ~Region();

    inline void setMarker(SIndex marker) { marker_ = marker; }

    inline SIndex marker() const { return marker_; }

    void resize(const Mesh & mesh, SIndex cellMarker);

    void resize(const std::vector < Cell * > & cells);

    inline void markBackground(bool background){ isBackground_ = background; }

    void setBackground(bool background = true);

    inline bool isBackground() const { return isBackground_ ; }

    inline void markSingle(bool issingle){ isSingle_ = issingle; }

    void setSingle(bool background = true);

    inline bool isSingle() const { return isSingle_ ; }

    inline bool isInParaDomain() const { return _isInParaDomain; }

    inline Index parameterCount() const { return parameterCount_; }

    inline Index startParameter() const { return startParameter_; }

    inline Index endParameter() const { return endParameter_; }

    void setConstraintType(Index type);

    inline Index constraintType() const { return constraintType_; }

    Index constraintCount() const;

    void fillConstraints(RSparseMapMatrix & C, Index startConstraintsID);

    void setConstraintWeights(double bc);

    void setConstraintWeights(const RVector & cw);

    const RVector & constraintWeights();

    void fillConstraintWeights(RVector & vec, Index constraintStart);

    void _createConstraintWeights();

    inline void setZWeight(double zw){
        zWeight_ = zw;
        this->constraintWeights_.clear();
    }

    inline double zWeight() const { return zWeight_; }

    void setFixValue(double val); 
    inline double fixValue() const { return fixValue_;}
 
    void fillBoundaryNorm(std::vector< RVector3 > & vnorm, Index boundCount);
 
    void fillBoundarySize(RVector & vec, Index boundStart);
 
    void fillStartModel(RVector & vec);
 
//     void fillStartVector(RVector & vec);
 
    void fillModelControl(RVector & vec);

    const std::vector < Boundary * > & boundaries() const { return bounds_; }
 
    const std::vector < Cell * > & cells() const { return cells_; }
 
    std::vector < Cell * > & cells() { return cells_; }
 
    void countParameter(Index start);

    bool isPermuted() const { return isPermuted_; }

    void permuteParameterMarker(const IndexArray & p);

    const IndexArray & paraIds() const { return paraIDs_; }

    void setStartModel(const RVector & start);

    void setStartModel(double start);
 
//     /*! DEPRECATED use setStartModel */
//     void setStartVector(const RVector & start);
//     /*! DEPRECATED use setStartModel */
    void setStartValue(double start){ DEPRECATED setStartModel(start);}
 
    void setModelControl(double val);
    void setModelControl(const RVector & vec);
 
    inline double modelControl() { return modelControl_; }
 
    void setTransModel(Trans< RVector > & tM);

    Trans< RVector > * transModel() { return tM_; }

    void setLowerBound(double lb);
    void setUpperBound(double ub);

    void setParameters(double start, double lb, double ub, std::string transString="");
 
    void setModelTransStr_(const std::string & val);
 
    void setModelControlStr_(  const std::string & val){ setModelControl(toDouble(val)); }
    void setStartModelStr_(    const std::string & val){ setStartModel(toDouble(val)); }
    void setZWeightStr_(       const std::string & val){ setZWeight(toDouble(val)); }
    void setFixValueStr_(      const std::string & val){ setFixValue(toDouble(val)); }
    void setConstraintTypeStr_(const std::string & val){ setConstraintType(toInt(val)); }
    void setLowerBoundStr_(    const std::string & val){ setLowerBound(toDouble(val)); }
    void setUpperBoundStr_(    const std::string & val){ setUpperBound(toDouble(val)); }
    void setSingleStr_(        const std::string & val){ markSingle(toInt(val) != 0); }
    void setBackgroundStr_(    const std::string & val){ markBackground(toInt(val) != 0); }
 
protected:
    void init_();
    void copy_(const Region & region);
 
    SIndex marker_;
    RegionManager * parent_;
 
    std::vector < Cell * > cells_;
    mutable std::vector < Boundary * > bounds_;
 
    bool isBackground_;
    bool isSingle_;
    bool isPermuted_;
    bool _isInParaDomain;
 
    IndexArray paraIDs_;
 
    Index parameterCount_;
    Index startParameter_;
    Index endParameter_;
 
    Index constraintType_;
 
    RVector startModel_;
    double modelControl_;
    RVector constraintWeights_;
 
    double zWeight_;
    double fixValue_;
 
    double mcDefault_;
    double startDefault_;
 
//    double minZWeight_;
 
    double lowerBound_;
    double upperBound_;
 
    Trans< RVector >        * tM_;
    bool ownsTrans_;        // smart ptr. would be nice here.
    std::string transString_;
};
 
class DLLEXPORT RegionManager{
public:
    RegionManager(bool verbose=true);
 
    RegionManager & operator = (const RegionManager & rm);
 
    ~RegionManager();
 
    void setVerbose(bool verbose){ verbose_ = verbose; }
    bool verbose() const { return verbose_; }
 
    void clear();
 
    void setMesh(const Mesh & mesh, bool holdRegionInfos=false);
 
    const Mesh & mesh() const;
 
    inline Mesh * pMesh() { return mesh_; }

    Region * addRegion(SIndex marker);
    
    Region * addRegion(SIndex marker, const Mesh & mesh){
        return addRegion(marker, mesh, marker);
    }
    Region * addRegion(SIndex marker, const Mesh & mesh, SIndex cellMarker);
#ifndef PYGIMLI_CAST
    const std::map < SIndex, Region * > & regions() const { return regionMap_; }
 
    std::map < SIndex, Region * >  * regions() { return &regionMap_; }
#endif
    Index regionCount() const { return regionMap_.size(); }

    IVector regionIdxs() const {
        return this->allRegionMarker_(false);
    }

    Region * region(SIndex marker);
 
    inline bool regionExists(SIndex marker) { return (regionMap_.count(marker) > 0); }

    void loadMap(const std::string & fname);

    void saveMap(const std::string & fname);

    inline void setParameterCount(Index count) { parameterCount_ = count; }

    Index parameterCount() const;

    Index constraintCount() const;
 
    Index interRegionConstraintsCount() const;

    RVector createStartModel();

    void fillStartModel(RVector & vec);

    void fillStartVector(RVector & vec){DEPRECATED;
         fillStartModel(vec);}

     RVector createStartVector(){DEPRECATED;
         return createStartModel();}

    RVector createModelControl();

    void fillModelControl(RVector & vec);

    RVector constraintWeights();

    void fillConstraintWeights(RVector & vec);

    void fillConstraints(RSparseMapMatrix & C);

    void setZWeight(double z);
 
    void setConstraintType(Index type);
 
    void fillBoundarySize(RVector & vec);
 
    inline const Mesh & paraDomain() const { return *paraDomain_; }
 
    inline Mesh & paraDomain() { return *paraDomain_; }
 
    std::vector < RVector3 > boundaryNorm() const;

    void permuteParameterMarker(const IVector & p);
 
    void createParaDomain_();
 
    void recountParaMarker_();
 
    TransCumulative < RVector > * transModel();
 
    void setLocalTransFlag(bool flag);
 
    bool haveLocalTrans() const { return haveLocalTrans_; }

    void setInterRegionConstraint(SIndex a, SIndex b, double weight);

    const std::map< std::pair< SIndex, SIndex >, double > interRegionConstraints() const { 
        return this->interRegionConstraints_; 
    }

    void setInterfaceConstraint(SIndex marker, double weight) { 
        this->interfaceConstraints_[marker] = weight; 
    }

    const std::map< SIndex, double > interfaceConstraints() const { 
        return this->interfaceConstraints_; 
    }
 
    void findInterRegionInterfaces();
    
protected:
    void copy_(const RegionManager & rm);


    Region * createRegion_(SIndex marker, const Mesh & mesh, SIndex cellMarker);

    Region * createSingleRegion_(SIndex marker, const std::vector < Cell * > & cells);
 
    IVector allRegionMarker_(bool excludeBoundary=false) const;
 
    bool verbose_;
    bool isPermuted_;
 
    Index parameterCount_;
 
    Mesh * mesh_;
    Mesh * paraDomain_;
 
    std::map < SIndex, Region * > regionMap_;
    std::map< std::pair< SIndex, SIndex >, std::list < Boundary * > > interRegionInterfaceMap_;
    std::map< std::pair< SIndex, SIndex >, double > interRegionConstraints_;
    std::map< SIndex, double > interfaceConstraints_;
 
    RVector _cWeights; // cache cWeights to avoid double creation.
 
    double interRegionConstraintZWeights_;
 
    TransCumulative < RVector > localTrans_;
    bool haveLocalTrans_;
    bool localTransHaveChanges_;
};
 
 
} // namespace GIMLI
 
#endif // _GIMLI_REGIONMANAGER__H