TPZMHMeshControl.h

Go to the documentation of this file.

//
//  TPZMHMeshControl.h
//  PZ
//
//  Created by Philippe Devloo on 5/3/14.
//
//

#ifndef __PZ__TPZMHMeshControl__
#define __PZ__TPZMHMeshControl__

#include <iostream>

#include "pzgmesh.h"
#include "pzcmesh.h"
#include "pzcompel.h"

class TPZMHMeshControl
{
    
public:

    enum MProblemType {ENone, EScalar, EElasticity2D, EElasticity3D};

protected:
    TPZAutoPointer<TPZGeoMesh> fGMesh;
    
    TPZAutoPointer<TPZCompMesh> fCMesh;
    
    TPZAutoPointer<TPZCompMesh> fCMeshLagrange;
    
    TPZAutoPointer<TPZCompMesh> fCMeshConstantPressure;
    
    // this mesh is the same as fCMesh if there are no lagrange multipliers assocated with the average pressure
    TPZAutoPointer<TPZCompMesh> fPressureFineMesh;
    
    MProblemType fProblemType = EScalar;
    
    int fNState = 1;
    
public:
    int fSkeletonMatId = 4;
    
    int fSecondSkeletonMatId = 3;
    
    int fPressureSkeletonMatId = 6;
    
    int fLagrangeMatIdLeft = 550, fLagrangeMatIdRight = 551;
    
    std::set<int> fMaterialIds;
    
    std::set<int> fMaterialBCIds;
    
protected:
    
    int fpOrderInternal = 1;
    
    int fpOrderSkeleton = 1;
    
    //internal order for enrichement spaces
    
    int fHdivmaismais = 0;
    
    int fSkeletonWrapMatId = 500;
    
    int fBoundaryWrapMatId = 502;
    
    int fInternalWrapMatId = 501;
    
    TPZManVector<int64_t> fGeoToMHMDomain;
    
    std::map<int64_t,int64_t> fMHMtoSubCMesh;
    
    // if fInterfaces[index].second == index, then the skeleton is boundary skeleton
    std::map<int64_t, std::pair<int64_t,int64_t> > fInterfaces;
    
    std::map<TPZCompMesh *,TPZManVector<int64_t> > fConnectToSubDomainIdentifier;
    

    bool fLagrangeAveragePressure = false;
    
    bool fHybridize = false;
    
    bool fSwitchLagrangeSign = false;
    
public:
    int64_t fGlobalSystemSize;
    
    int64_t fGlobalSystemWithLocalCondensationSize;
    
    int64_t fNumeq;
public:
    TPZMHMeshControl() : fProblemType(EScalar), fNState(1), fGeoToMHMDomain(), fMHMtoSubCMesh(),
    fLagrangeAveragePressure(false),
    fHybridize(false), fSwitchLagrangeSign(false), fGlobalSystemSize(-1), fGlobalSystemWithLocalCondensationSize(-1), fNumeq(-1)
    {
        
    }

    TPZMHMeshControl(int dimension) : fProblemType(EScalar), fNState(1), fGeoToMHMDomain(), fMHMtoSubCMesh(),
    fLagrangeAveragePressure(false),
    fHybridize(false), fSwitchLagrangeSign(false), fGlobalSystemSize(-1), fGlobalSystemWithLocalCondensationSize(-1), fNumeq(-1)
    {
        fGMesh = new TPZGeoMesh;
        fGMesh->SetDimension(dimension);
        fCMesh = new TPZCompMesh(fGMesh);
        fPressureFineMesh = new TPZCompMesh(fGMesh);
    }

    //    TPZMHMeshControl(TPZAutoPointer<TPZGeoMesh> gmesh, std::set<int64_t> &coarseindices);
    
    TPZMHMeshControl(TPZAutoPointer<TPZGeoMesh> gmesh, TPZVec<int64_t> &geotomhm);
    
    TPZMHMeshControl(TPZAutoPointer<TPZGeoMesh> gmesh);
    
    TPZMHMeshControl(const TPZMHMeshControl &copy);
    
    virtual ~TPZMHMeshControl();
    
    TPZMHMeshControl &operator=(const TPZMHMeshControl &cp);
    
    TPZAutoPointer<TPZCompMesh> CMesh() const
    {
        return fCMesh;
    }
    
    TPZAutoPointer<TPZGeoMesh> GMesh() const
    {
        return fGMesh;
    }
    
    TPZAutoPointer<TPZCompMesh> PressureMesh()
    {
        return fPressureFineMesh;
    }
    
    void DefinePartitionbyCoarseIndices(TPZVec<int64_t> &coarseindices);

    void DefineSkeleton(std::map<int64_t,std::pair<int64_t,int64_t> > &skeleton)
    {
        fInterfaces = skeleton;
#ifdef PZDEBUG
        // the skeleton elements link MHM domain indices (which need to exist)
        for (auto it = fInterfaces.begin(); it != fInterfaces.end(); it++) {
            if (fMHMtoSubCMesh.find(it->second.first) == fMHMtoSubCMesh.end()) {
                DebugStop();
            }
            if (it->second.second != it->first && fMHMtoSubCMesh.find(it->second.second) == fMHMtoSubCMesh.end()) {
                DebugStop();
            }
        }
#endif
    }
    
    void DefinePartition(TPZVec<int64_t> &partitionindex, std::map<int64_t,std::pair<int64_t,int64_t> > &skeleton);

    // This method calculates the skeleton indexes
    void DefinePartition(TPZVec<int64_t> &partitionindex);

    void SetProblemType(MProblemType problem)
    {
        fProblemType = problem;
        switch(problem)
        {
            case EScalar:
                fNState = 1;
                break;
            case EElasticity2D:
                fNState = 2;
                break;
            case EElasticity3D:
                fNState = 3;
                break;
            default:
                DebugStop();
        }
    }
    void SetInternalPOrder(int order)
    {
        fpOrderInternal = order;
        if(fCMesh) fCMesh->SetDefaultOrder(order);
    }
    
    void SetSkeletonPOrder(int order)
    {
        fpOrderSkeleton = order;
    }
    
    
    void SetHdivmaismaisPOrder(int order)
    {
        fHdivmaismais = order;
    }

    
    void SetLagrangeAveragePressure(bool flag)
    {
        fLagrangeAveragePressure = flag;
    }
    
    virtual void SetHybridize(bool flag = true)
    {
        fHybridize = flag;
    }
    
    virtual bool GetHybridize() {
        return fHybridize;
    }
    virtual void BuildComputationalMesh(bool usersubstructure);
    
protected:
    virtual void CreateSkeletonElements();
    
    virtual void InsertPeriferalMaterialObjects();
public:
    void DivideSkeletonElements(int ndivide);
    
    void DivideBoundarySkeletonElements();
    
    void SwitchLagrangeMultiplierSign(bool sw)
    {
        fSwitchLagrangeSign = sw;
    }
    
    void Print(std::ostream &out);
    
    void PrintDiagnostics(std::ostream &out);
    
    void GetMeshVec(TPZVec<TPZCompMesh *> &meshvec)
    {
        if (fCMeshLagrange)
        {
            meshvec.Resize(3);
            meshvec[0] = fPressureFineMesh.operator->();
            meshvec[1] = fCMeshLagrange.operator->();
            meshvec[2] = fCMeshConstantPressure.operator->();
        }
        else
        {
            meshvec.resize(0);
        }
    }
    
    TPZVec<TPZAutoPointer<TPZCompMesh> > GetMeshes()
    {
        TPZManVector<TPZAutoPointer<TPZCompMesh>,3> result;
        if (fCMeshLagrange)
        {
            result.Resize(3);
            result[0] = fPressureFineMesh;
            result[1] = fCMeshLagrange;
            result[2] = fCMeshConstantPressure;
        }
        else
        {
            result.resize(0);
        }
        return result;
    }

    TPZManVector<int64_t> GetGeoToMHMDomain() {
        return fGeoToMHMDomain;
    }

    std::map<int64_t,int64_t> &Coarse_to_Submesh()
    {
        return fMHMtoSubCMesh;
    }
    
    //  only implemented in the hybrid Hdiv version
    virtual void CheckMeshConsistency()
    {
        
    }
    
protected:
    TPZCompMesh *CriaMalhaTemporaria();
    
    void CreateInternalElements();
    
    void AddBoundaryElements();
    
    void AddBoundaryInterfaceElements();
    
    void CreateSkeleton();
    
    void CreateInterfaceElements();
    void CreateInterfaceElements2();
    
    virtual void HybridizeSkeleton(int skeletonmatid, int pressurematid);
    
    bool IsSibling(int64_t son, int64_t father);
    
    void AddElementBoundaries(int64_t elseed, int64_t compelindex, TPZStack<TPZCompElSide> &result);
    
    void TransferToMultiphysics();
    
    void SubStructure();
    void SubStructure2();
    
    void BuildWrapMesh(int dim);
    
    int HasWrapNeighbour(TPZGeoElSide gelside);
    
    void CheckDivisionConsistency(TPZGeoElSide gelside);
    
    int WrapMaterialId(TPZGeoElSide gelside);
    
    virtual bool IsSkeletonMatid(int matid)
    {
        return matid == fSkeletonMatId;
    }
    
    virtual bool IsBoundaryMatid(int matid)
    {
        return fMaterialBCIds.find(matid) != fMaterialBCIds.end();
    }
    
    void CreateWrap(TPZGeoElSide gelside);
    
    void CreateWrap(TPZGeoElSide gelside, int wrapmaterial);
    
    void DivideWrap(TPZGeoEl *wrapelement);
    
protected:
    void CreateLagrangeMultiplierMesh();
    
    void SetSubdomain(TPZCompEl *cel, int64_t subdomain);
    
    void SetSubdomain(TPZCompMesh *cmesh, int64_t connectindex, int64_t subdomain);
    
    // this method calls debugstop if the element beint64_ts to two subdomains
    int64_t WhichSubdomain(TPZCompEl *cel);
    
    // the subdomain indices of the connects to the multi physics mesh
    void JoinSubdomains(TPZVec<TPZCompMesh *> &meshvec, TPZCompMesh *multiphysicsmesh);
    
    void PrintSubdomain(int64_t elindex, std::ostream &out);
    
    void PrintBoundaryInfo(std::ostream &out);
    
    // the computational mesh is determined by the element pointed to by the geometric element
    void ConnectedElements(int64_t skeleton, std::pair<int64_t,int64_t> &leftright, std::map<int64_t, std::list<TPZCompElSide> > &ellist);

    // the computational mesh is determined by the element pointed to by the geometric element
    void ConnectedInterfaceElements(int64_t skeleton, std::pair<int64_t,int64_t> &leftright, std::map<int64_t, std::list<TPZInterfaceElement *> > &ellist);
};

#endif /* defined(__PZ__TPZMHMeshControl__) */