TPZSSpStructMatrix.cpp

Go to the documentation of this file.

#include "TPZSSpStructMatrix.h"
#include "pzstrmatrix.h"

#include "pzgeoelbc.h"
#include "pzgmesh.h"
#include "pzcmesh.h"

#include "pzanalysis.h"
#include "pzsolve.h"
#include "pzstepsolver.h"

#include "pzdxmesh.h"
#include <fstream>

#include "pzelmat.h"

#include "pzsysmp.h"
#include "pzmetis.h"
#include "pzbndcond.h"
#include "TPZTimer.h"

#include "pzlog.h"
#ifdef LOG4CXX
static LoggerPtr logger(Logger::getLogger("pz.StrMatrix"));
#endif

using namespace std;

TPZStructMatrix * TPZSymetricSpStructMatrix::Clone(){
    return new TPZSymetricSpStructMatrix(*this);
}
TPZMatrix<STATE> * TPZSymetricSpStructMatrix::CreateAssemble(TPZFMatrix<STATE> &rhs,
                                              TPZAutoPointer<TPZGuiInterface> guiInterface){
                
#ifdef LOG4CXX
    if (logger->isDebugEnabled())
    {
        LOGPZ_DEBUG(logger,"TPZSymetricSpStructMatrix::CreateAssemble starting");
    }
#endif
                
    int64_t neq = fMesh->NEquations();
    if(fMesh->FatherMesh()) {
                                cout << "TPZSymetricSpStructMatrix should not be called with CreateAssemble for a substructure mesh\n";
                                return new TPZSYsmpMatrix<STATE>(0,0);
    }
//    std::cout << "Creating\n";
    TPZMatrix<STATE> *stiff = Create();//new TPZFYsmpMatrix(neq,neq);
    TPZSYsmpMatrix<STATE> *mat = dynamic_cast<TPZSYsmpMatrix<STATE> *> (stiff);
    rhs.Redim(neq,1);
    //stiff->Print("Stiffness TPZFYsmpMatrix :: CreateAssemble()");
    TPZTimer before("Assembly of a sparse matrix");
//    std::cout << "Assembling\n";
    before.start();
#ifdef LOG4CXX
    if(logger->isDebugEnabled()) LOGPZ_DEBUG(logger,"TPZSymetricSpStructMatrix::CreateAssemble calling Assemble()");
#endif
                Assemble(*stiff,rhs,guiInterface);
    mat->ComputeDiagonal();
    
//    std::cout << "Rhs norm " << Norm(rhs) << std::endl;
    
    before.stop();
    //std::cout << __PRETTY_FUNCTION__ << " " << before << std::endl;
    //    mat->ComputeDiagonal();
    //stiff->Print("Stiffness TPZFYsmpMatrix :: CreateAssemble()");
#ifdef LOG4CXX
    if(logger->isDebugEnabled()) LOGPZ_DEBUG(logger,"TPZSymetricSpStructMatrix::CreateAssemble exiting");
#endif
    return stiff;
}
TPZMatrix<STATE> * TPZSymetricSpStructMatrix::Create(){

    TPZStack<int64_t> elgraph;
    TPZVec<int64_t> elgraphindex;
    //    int nnodes = 0;
    fMesh->ComputeElGraph(elgraph,elgraphindex);
    TPZMatrix<STATE> * mat = SetupMatrixData(elgraph, elgraphindex);
    return mat;
}

TPZMatrix<STATE> * TPZSymetricSpStructMatrix::SetupMatrixData(TPZStack<int64_t> & elgraph, TPZVec<int64_t> &elgraphindex){
    
    int64_t neq = fEquationFilter.NActiveEquations();
    TPZSYsmpMatrix<STATE> * mat = new TPZSYsmpMatrix<STATE>(neq,neq);
    
    TPZMetis metis;
    metis.SetElementsNodes(elgraphindex.NElements() -1 ,fMesh->NIndependentConnects());
    metis.SetElementGraph(elgraph,elgraphindex);
    
    TPZManVector<int64_t> nodegraph;
    TPZManVector<int64_t> nodegraphindex;
    metis.ConvertGraph(elgraph,elgraphindex,nodegraph,nodegraphindex);
    int64_t i;
    int64_t nblock = nodegraphindex.NElements()-1;
    // number of values in the sparse matrix
    int64_t totalvar = 0;
    // number of equations
    int64_t totaleq = 0;
    for(i=0;i<nblock;i++){
        int64_t iblsize = fMesh->Block().Size(i);
        int64_t iblpos = fMesh->Block().Position(i);
        int64_t numactive = fEquationFilter.NumActive(iblpos, iblpos+iblsize);
        if (!numactive) {
            continue;
        }
        totaleq += iblsize;
        int64_t icfirst = nodegraphindex[i];
        int64_t iclast = nodegraphindex[i+1];
        int64_t j;
        //longhin
        totalvar+=(iblsize*(iblsize+1))/2;
        for(j=icfirst;j<iclast;j++) {
            int64_t col = nodegraph[j];
            if (col < i) {
                continue;
            }
            
            if (col == i) {
                DebugStop();
            }
            
            int64_t colsize = fMesh->Block().Size(col);
            int64_t colpos = fMesh->Block().Position(col);
            int64_t numactive = fEquationFilter.NumActive(colpos, colpos+colsize);
            if (!numactive) {
                continue;
            }
            totalvar += iblsize*colsize;
        }
    }
    
    int64_t ieq = 0;
    // pos is the position where we will put the column value
    int64_t pos = 0;
    
    nblock=fMesh->NIndependentConnects();
    
    TPZVec<int64_t> Eq(totaleq+1);
    TPZVec<int64_t> EqCol(totalvar);
    TPZVec<STATE> EqValue(totalvar,0.);
    for(i=0;i<nblock;i++){
        int64_t iblsize = fMesh->Block().Size(i);
        int64_t iblpos = fMesh->Block().Position(i);
        TPZManVector<int64_t> rowdestindices(iblsize);
        for (int64_t i=0; i<iblsize; i++) {
            rowdestindices[i] = iblpos+i;
        }
        fEquationFilter.Filter(rowdestindices);

        int64_t ibleq;
        // working equation by equation
        for(ibleq=0; ibleq<rowdestindices.size(); ibleq++) {
            if (rowdestindices[ibleq] != ieq) {
                DebugStop();
            }
            Eq[ieq] = pos;
            int64_t colsize,colpos,jbleq;
            int64_t diagonalinsert = 0;
            int64_t icfirst = nodegraphindex[i];
            int64_t iclast = nodegraphindex[i+1];
            int64_t j;
            for(j=icfirst;j<iclast;j++)
            {
                int64_t col = nodegraph[j];
                if (col < i) {
                    continue;
                }
                // force the diagonal block to be inserted
                // the nodegraph does not contain the pointer to itself
                if(!diagonalinsert && col > i)
                {
                    diagonalinsert = 1;
                    int64_t colsize = fMesh->Block().Size(i);
                    int64_t colpos = fMesh->Block().Position(i);
                    TPZManVector<int64_t> destindices(colsize);
                    for (int64_t i=0; i<colsize; i++) {
                        destindices[i] = colpos+i;
                    }
                    fEquationFilter.Filter(destindices);
                    int64_t jbleq;
                    for(jbleq=0; jbleq<destindices.size(); jbleq++) {
                        //             if(colpos+jbleq == ieq) continue;
                        int64_t jeq = destindices[jbleq];
                        if (jeq < ieq) {
                            continue;
                        }
                        EqCol[pos] = destindices[jbleq];
                        EqValue[pos] = 0.;
                        //            colpos++;
                        pos++;
                    }
                }
                colsize = fMesh->Block().Size(col);
                colpos = fMesh->Block().Position(col);
                if (fEquationFilter.NumActive(colpos, colpos+colsize) == 0) {
                    continue;
                }
                TPZManVector<int64_t> destindices(colsize);
                for (int64_t i=0; i<colsize; i++) {
                    destindices[i] = colpos+i;
                }
                fEquationFilter.Filter(destindices);
                for(jbleq=0; jbleq<destindices.size(); jbleq++) {
                    int64_t jeq = destindices[jbleq];
                    if (jeq < ieq) {
                        continue;
                    }
                    EqCol[pos] = jeq;
                    EqValue[pos] = 0.;
                    colpos++;
                    pos++;
                }
            }
            // all elements are below (last block certainly)
            if(!diagonalinsert)
            {
                diagonalinsert = 1;
                int64_t colsize = fMesh->Block().Size(i);
                int64_t colpos = fMesh->Block().Position(i);
                TPZManVector<int64_t> destindices(colsize);
                for (int64_t i=0; i<colsize; i++) {
                    destindices[i] = colpos+i;
                }
                fEquationFilter.Filter(destindices);
                int64_t jbleq;
                for(jbleq=0; jbleq<destindices.size(); jbleq++) {
                    //             if(colpos+jbleq == ieq) continue;
                    int64_t jeq = destindices[jbleq];
                    if (jeq < ieq) {
                        continue;
                    }
                    EqCol[pos] = jeq;
                    EqValue[pos] = 0.;
                    //            colpos++;
                    pos++;
                }
            }
            ieq++;
        }
    }

    Eq[ieq] = pos;
    mat->SetData(Eq,EqCol,EqValue);
    return mat;
}

TPZSymetricSpStructMatrix::TPZSymetricSpStructMatrix() : TPZStructMatrix(){
    
}

TPZSymetricSpStructMatrix::TPZSymetricSpStructMatrix(TPZCompMesh *mesh) : TPZStructMatrix(mesh)
{}

#ifndef STATE_COMPLEX
#include "pzmat2dlin.h"

int TPZSymetricSpStructMatrix::main() {
                
                int refine=5;
                int order=5;
                
                TPZGeoMesh gmesh;
                TPZCompMesh cmesh(&gmesh);
                double coordstore[4][3] = {{0.,0.,0.},{1.,0.,0.},{1.,1.,0.},
                                {0.,1.,0.}};
                
                int i,j;
                TPZVec<REAL> coord(3,0.);
                for(i=0; i<4; i++) {
                                // initializar as coordenadas do no em um vetor
                                for (j=0; j<3; j++) coord[j] = coordstore[i][j];
                                
                                // identificar um espa� no vetor onde podemos armazenar
                                // este vetor
                                gmesh.NodeVec ().AllocateNewElement ();
                                
                                // initializar os dados do n�       gmesh.NodeVec ()[nodeindex].Initialize (i,coord,gmesh);
                }
                int el;
                TPZGeoEl *gel;
                for(el=0; el<1; el++) {
                                
                                // initializar os indices dos nos
                                TPZVec<int64_t> indices(4);
                                for(i=0; i<4; i++) indices[i] = i;
                                // O proprio construtor vai inserir o elemento na malha
                                //       gel = new TPZGeoElQ2d(el,indices,1,gmesh);
                                int64_t index;
                                gel = gmesh.CreateGeoElement(EQuadrilateral,indices,1,index);
                }
                gmesh.BuildConnectivity ();
                
                TPZVec<TPZGeoEl *> subel;
                
                cout << "Refinement ";
                cin >> refine;
                cout << endl;
                DebugStop();
//              UniformRefine(refine,gmesh);
                
                TPZGeoElBC gelbc(gel,4,-4);
                TPZMat2dLin *meumat = new TPZMat2dLin(1);
                TPZFMatrix<STATE> xk(1,1,1.),xc(1,2,0.),xf(1,1,1.);
                meumat->SetMaterial (xk,xc,xf);
                TPZMaterial * meumatptr(meumat);
                cmesh.InsertMaterialObject(meumatptr);
                
                TPZFMatrix<STATE> val1(1,1,0.),val2(1,1,0.);
                TPZMaterial * bnd = meumat->CreateBC (meumatptr,-4,0,val1,val2);
                cmesh.InsertMaterialObject(bnd);
                
                cout << "Interpolation order ";
                cin >> order;
                cout << endl;
                
                //TPZCompEl::gOrder = order;
                cmesh.SetDefaultOrder(order);
                
                cmesh.AutoBuild();
                //              cmesh.AdjustBoundaryElements();
                cmesh.InitializeBlock();
                
                ofstream output("outputPar.dat");
                TPZAnalysis an(&cmesh,true,output);
                
                //TPZVec<int> numelconnected(cmesh.NEquations(),0);
                TPZSymetricSpStructMatrix mat(&cmesh);
                
                an.SetStructuralMatrix(mat);
                
                TPZStepSolver<STATE> sol;
                sol.SetJacobi(100,1.e-5,0);
                
                
                an.SetSolver(sol);
                an.Run(output);
                output.flush();
                cout.flush();
                return 0;
                
}
#endif