dcfemmodelling.h

/******************************************************************************
 *   Copyright (C) 2006-2024 by the resistivity.net development team          *
 *   Carsten Rücker carsten@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 _BERT_DCFEMMODDELING__H
#define _BERT_DCFEMMODDELING__H
 
#include "bert.h"
#include "datamap.h"
 
#include <modellingbase.h>
#include <sparsematrix.h>
#include <pos.h>
#include <matrix.h>
 
#include <vector>
 
namespace GIMLI{

DLLEXPORT void dcfemDomainAssembleStiffnessMatrix(RSparseMatrix & S, const Mesh & mesh,
                                                  double k=0.0, bool fix=true);
 
DLLEXPORT void dcfemDomainAssembleStiffnessMatrix(CSparseMatrix & S, const Mesh & mesh,
                                                  double k=0.0, bool fix=true);
 
// DLLEXPORT void assembleStiffnessMatrixHomogenDirichletBC(RSparseMatrix & S,
//                                                          const IndexArray & nodeID);
 
DLLEXPORT double mixedBoundaryCondition(const Boundary & boundary,
                                        const RVector3 & sourcePos,
                                        double k=0.0);
 
DLLEXPORT void assembleCompleteElectrodeModel(RSparseMatrix & S,
                                              const std::vector < ElectrodeShape * > & elecs,                                               
                                              uint oldMatSize, bool lastIsReferenz,
                                              const RVector & contactImpedances);

DLLEXPORT void setComplexResistivities(Mesh & mesh, const CVector & res);

DLLEXPORT void setComplexResistivities(Mesh & mesh,
                                      const std::map < float, Complex > & aMap);

DLLEXPORT void setComplexResistivities(Mesh & mesh,
                                       const RVector & amp,
                                       const RVector & phase);

DLLEXPORT CVector getComplexResistivities(const Mesh & mesh);

DLLEXPORT void setComplexData(DataContainer & data, const CVector & u);

DLLEXPORT void setComplexData(DataContainer & data,
                              const RVector & re, const RVector & im);

DLLEXPORT CVector getComplexData(const DataContainer & data);
 
 
class DLLEXPORT DCMultiElectrodeModelling : public GIMLI::ModellingBase {
public:
    DCMultiElectrodeModelling(bool verbose=false);
 
    DCMultiElectrodeModelling(Mesh & mesh,
                              bool verbose=false);
 
    DCMultiElectrodeModelling(DataContainerERT & dataContainer,
                              bool verbose=false);
 
    DCMultiElectrodeModelling(Mesh & mesh, DataContainerERT & dataContainer,
                              bool verbose=false);
 
    virtual ~DCMultiElectrodeModelling();

 
    void assembleStiffnessMatrixDCFEMByPass(RSparseMatrix & S);
 
    void assembleStiffnessMatrixDCFEMByPass(CSparseMatrix & S);
 
    virtual void createJacobian(const RVector & model);
 
    virtual void createConstraints();
 
    DataContainerERT & dataContainer() const ;

    virtual RVector createDefaultStartModel();
 
//     RVector operator () (const RVector & model, double background=0) {
//         return response(model, background); }

    void mapERTModel(const CVector & model, Complex background);

    void mapERTModel(const RVector & model, double background);

    virtual RVector response(const RVector & model){return response(model, -9e99);}

    RVector response(const RVector & model, double background);
 
    void createCurrentPattern(std::vector < ElectrodeShape * > & eA,
                              std::vector < ElectrodeShape * > & eB,
                              bool reciprocity);
 
    virtual void calculate(DataMap & dMap);
 
    virtual void calculate(DataContainerERT & data, bool reciprocity=false);
 
    virtual void calculate(const std::vector < ElectrodeShape * > & eA,
                           const std::vector < ElectrodeShape * > & eB);
 
    virtual void calculateK(const std::vector < ElectrodeShape * > & eA,
                            const std::vector < ElectrodeShape * > & eB,
                            RMatrix & solutionK, int kIdx);
 
    virtual void calculateK(const std::vector < ElectrodeShape * > & eA,
                            const std::vector < ElectrodeShape * > & eB,
                            CMatrix & solutionK, int kIdx);
 
    template < class ValueType >
    void calculateKAnalyt(const std::vector < ElectrodeShape * > & eA,
                          const std::vector < ElectrodeShape * > & eB,
                          Matrix < ValueType > & solutionK,
                          double k, int kIdx) const;
 
    void setTopography(bool topography) { topography_=topography; }
    bool topography() const { return topography_; }
 
    bool neumann() const { return neumannDomain_; }

    void setComplex(bool c);

    bool complex() const { return complex_; }

    void setDipoleCurrentPattern(bool dipole) {dipoleCurrentPattern_=dipole;}

    bool dipoleCurrentPattern() const { return dipoleCurrentPattern_; }

    inline void setAnalytical(bool ana){ analytical_=ana; }
    inline bool analytical() const { return analytical_; }
 
    void collectSubPotentials(RMatrix & subSolutions){
        subSolutions_=& subSolutions;
    }

    inline const std::map < Index, Index > & currentPatternIdxMap() const {
        return currentPatternIdxMap_;
    }

    inline void setBypassMapFile(const std::string & fileName){
        byPassFile_=fileName;
    }
 
    inline void setkValues(const RVector & v) { kValues_=v; }
    inline const RVector & kValues() const { return kValues_; }
 
    inline void setWeights(const RVector & v) { weights_=v; }
    inline const RVector & weights() const { return weights_; }
 
    inline const std::vector< ElectrodeShape * > & electrodes() const {
        return electrodes_;
    }
 
    void setContactImpedances(const RVector & zi);
 
    virtual RVector calcGeometricFactor(const DataContainerERT & data,
                                        Index nModel=0);
 
    virtual void preCalculate(const std::vector < ElectrodeShape * > & eA,
                              const std::vector < ElectrodeShape * > & eB){
    }

    void setSingValue(bool s=true) { setSingValue_=s;}

    bool isSetSingValue() const { return setSingValue_;}

    void setSolver(SolverWrapper *solver){ solver_ = solver; }
    
private:
    void init_();
 
protected:
 
    template < class ValueType >
    void calculateK_(const std::vector < ElectrodeShape * > & eA,
                     const std::vector < ElectrodeShape * > & eB,
                     Matrix < ValueType > & solutionK, int kIdx);
 
    template < class ValueType >
    void assembleStiffnessMatrixDCFEMByPass_(SparseMatrix < ValueType > & S);
 
    template < class ValueType >
    DataMap response_(const Vector < ValueType > & model,
                                   ValueType background);
 
    template < class ValueType >
    Matrix < ValueType > * prepareJacobianT_(const Vector< ValueType > & model);
 
    RMatrix * prepareJacobian_(const RVector & model);
    CMatrix * prepareJacobian_(const CVector & model);
 
    void createJacobian_(const RVector & model, const RMatrix & u, RMatrix * J);
    void createJacobian_(const CVector & model, const CMatrix & u, CMatrix * J);
 
    virtual void deleteMeshDependency_();
    virtual void updateMeshDependency_();
    virtual void updateDataDependency_();
 
    virtual void searchElectrodes_();
 
    MatrixBase * subSolutions_;
 
    bool complex_;
 
    bool JIsRMatrix_;
    bool JIsCMatrix_;
 
    bool analytical_;
    bool topography_;
    bool neumannDomain_;
    bool subpotOwner_;
    bool lastIsReferenz_;
    bool setSingValue_;
 
    std::string byPassFile_;
 
    RVector kValues_;
    RVector weights_;
 
    double surfaceZ_;
 
    IndexArray calibrationSourceIdx_;
    IndexArray bypassNodeIdx_;
 
    std::vector< ElectrodeShape * > electrodes_;
    ElectrodeShape * electrodeRef_;
 
    std::vector< ElectrodeShape * > passiveCEM_;
 
    RVector3 sourceCenterPos_;
 
    bool buildCompleteElectrodeModel_;
 
    bool dipoleCurrentPattern_;
    std::map < Index, Index > currentPatternIdxMap_;
 
    RMatrix potentialsCEM_;
    RVector vContactImpedance_;
 
    DataMap * primDataMap_;
 
    SolverWrapper *solver_;
};
 
class DLLEXPORT DCSRMultiElectrodeModelling : public DCMultiElectrodeModelling {
public:
    DCSRMultiElectrodeModelling(bool verbose=false)
        : DCMultiElectrodeModelling(verbose) {
        init_();
    }
 
    DCSRMultiElectrodeModelling(Mesh & mesh, bool verbose=false)
        : DCMultiElectrodeModelling(mesh, verbose) {
        init_();
    }
 
    DCSRMultiElectrodeModelling(Mesh & mesh, DataContainerERT & dataContainer, bool verbose=false)
        : DCMultiElectrodeModelling(mesh, dataContainer, verbose){
        init_();
    }
 
    virtual ~DCSRMultiElectrodeModelling(){
        if (primPot_  && primPotOwner_)  delete primPot_;
        if (primMesh_ && primMeshOwner_) delete primMesh_;
    };
 
    virtual void calculateK(const std::vector < ElectrodeShape * > & eA,
                            const std::vector < ElectrodeShape * > & eB,
                            RMatrix & solutionK, int kIdx);
 
    inline void setPrimaryPotFileBody(const std::string & primPotFileBody){
        primPotFileBody_=primPotFileBody;
    }
 
    inline void setPrimaryPotential(RMatrix & primPot) { primPot_=& primPot; }
 
    inline RMatrix & primaryPotential() { return *primPot_; }
 
    inline void setPrimaryMesh(Mesh & mesh){ primMesh_=&mesh; }
 
    void setPrimaryMesh(const std::string & meshname);
 
    inline Mesh & primaryMesh() {return *primMesh_; }
 
    //const DataMap & primDataMap() const { return ; }
 
    virtual void preCalculate(const std::vector < ElectrodeShape * > & eA,
                               const std::vector < ElectrodeShape * > & eB);
 
private:
    void init_(){
        primPotFileBody_= NOT_DEFINED;
 
        primPotOwner_  =false;
        primPot_       =NULL;
 
        primMeshOwner_ =false;
        primMesh_      =NULL;
    }
 
protected:
    virtual void updateMeshDependency_();
    virtual void updateDataDependency_();
 
    void checkPrimpotentials_(const std::vector < ElectrodeShape * > & eA,
                               const std::vector < ElectrodeShape * > & eB);
 
    std::string primPotFileBody_;
 
    bool primPotOwner_;
    RMatrix * primPot_;
 
    bool primMeshOwner_;
    Mesh * primMesh_;
    Mesh mesh1_;
};
 
} //namespace BERT
 
#endif // _BERT_DCFEMMODDELING__H