substruct.cpp

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <pz_config.h>
#endif

#include <iostream>
#include <cstdlib>

#include "tpzdohrsubstruct.h"
#include "tpzdohrmatrix.h"
#include "tpzdohrprecond.h"
#include "pzdohrstructmatrix.h"
#include "pzstepsolver.h"
#include "pzcompel.h"
#include "pzgeoelbc.h"

#include "pzelast3d.h"
#include "pzbndcond.h"

#include "tpzdohrassembly.h"

#include "pzlog.h"
#include "tpzgensubstruct.h"
#include "tpzpairstructmatrix.h"
#include "pzviscoelastic.h"
#include "TPZVTKGeoMesh.h"
#include "pzgeotetrahedra.h"
#include "pzskylstrmatrix.h"

#include "tpzarc3d.h"
#include "tpzdohrmatrix.h"

#include "pzvtkmesh.h"

#include "pzlog.h"

#include "pzbfilestream.h" // TPZBFileStream, TPZFileStream
#include "pzmd5stream.h"

#include <fstream>
#include <string>
#include "TPZfTime.h"


#ifdef LOG4CXX
static LoggerPtr loggerconverge(Logger::getLogger("pz.converge"));
static LoggerPtr logger(Logger::getLogger("main"));
#endif

//#include "timing_analysis.h"
#include "arglib.h"
#include "run_stats_table.h"

#ifdef HAS_GETRUSAGE
#include <sys/resource.h> // getrusage
#endif

#ifdef USING_TBB
#include "tbb/task_scheduler_init.h"
using namespace tbb;
// If you have issues with: dyld: Library not loaded: libtbb.dylib
// try setting the LD path. Ex:
//   export DYLD_FALLBACK_LIBRARY_PATH=/Users/borin/Desktop/neopz/tbb40_297oss/lib/
#endif

void InsertElasticity(TPZAutoPointer<TPZCompMesh> mesh);
void InsertViscoElasticity(TPZAutoPointer<TPZCompMesh> mesh);
void InsertViscoElasticityCubo(TPZAutoPointer<TPZCompMesh> mesh);
TPZGeoMesh *MalhaPredio();
TPZGeoMesh *MalhaCubo();
void SetPointBC(TPZGeoMesh *gr, TPZVec<REAL> &x, int bc);
REAL Height(TPZGeoMesh *gmesh);
int SubStructure(TPZAutoPointer<TPZCompMesh> cmesh, REAL height);

using namespace std;

void help(const char* prg)
{
    cout << "Compute ...." << endl;
    cout << "The application is divided in four main steps: s1, s2, s3 and s4" << endl;
    cout << endl;
    cout << "  Step 1 (mp/mc  -> ckpt 1): " << endl;
    cout << "  Step 2 (ckpt 1 -> ckpt 2): dohrstruct->Create()" << endl;
    cout << "  Step 3 (ckpt 2 -> ckpt 3): dohrstruct->Assemble(), dohrstruct->Preconditioner() " << endl;
    cout << "  Step 4 (ckpt 3 -> end   ): cg.solve() ..." << endl << endl;
    cout << "Usage: " << prg << "starting_point stop_point file [-of output_file]" << endl << endl;
    cout << " starting_point = {-cf1|-cf2|-cf3|-mc|-mp} input_file" << endl;
    cout << " starting_point = {-st1|-st2|-st3}" << endl;
    
    clarg::arguments_descriptions(cout, "  ", "\n");
}

clarg::argString cf1("-cf1", "starts execution from checkpoint 1 (read checkpoint file)", "ckpt1.ckpt");
clarg::argString cf2("-cf2", "starts execution from checkpoint 2 (read checkpoint file)", "ckpt2.ckpt");
clarg::argString cf3("-cf3", "starts execution from checkpoint 3 (read checkpoint file)", "ckpt3.ckpt");
clarg::argBool   st1("-st1", "stop at checkpoint 1 (after dump)", false);
clarg::argBool   st2("-st2", "stop at checkpoint 2 (after dump)", false);
clarg::argBool   st3("-st3", "stop at checkpoint 3 (after dump)", false);

clarg::argString gen_sig_ckpt1("-gen_c1_md5", "generates MD5 signature for checkpoint 1 and dump into file.", "ckpt1.md5");
clarg::argString chk_sig_ckpt1("-chk_c1_md5", "compute MD5 signature for checkpoint 1 and check against MD5 at file.", "ckpt1.md5");
clarg::argString gen_sig_ckpt2("-gen_c2_md5", "generates MD5 signature for checkpoint 2 and dump into file.", "ckpt2.md5");
clarg::argString chk_sig_ckpt2("-chk_c2_md5", "compute MD5 signature for checkpoint 2 and check against MD5 at file.", "ckpt2.md5");
clarg::argString gen_sig_ckpt3("-gen_c3_md5", "generates MD5 signature for checkpoint 3 and dump into file.", "ckpt3.md5");
clarg::argString chk_sig_ckpt3("-chk_c3_md5", "compute MD5 signature for checkpoint 3 and check against MD5 at file.", "ckpt3.md5");
clarg::argString gen_sig_ckpt4("-gen_c4_md5", "generates MD5 signature for checkpoint 4 and dump into file.", "ckpt4.md5");
clarg::argString chk_sig_ckpt4("-chk_c4_md5", "compute MD5 signature for checkpoint 4 and check against MD5 at file.", "ckpt4.md5");

clarg::argString dc1("-dc1", "dump checkpoint 1 to file", "ckpt1.ckpt");
clarg::argString dc2("-dc2", "dump checkpoint 2 to file", "ckpt2.ckpt");
clarg::argString dc3("-dc3", "dump checkpoint 3 to file", "ckpt3.ckpt");
clarg::argString mc("-mc", "starts execution from beginning - read a \"malha_cubo\" input file",
                    "../cube1.txt");
clarg::argString mp("-mp", "starts execution from beginning - read a \"malha_predio\" input file",
                    "../8andares02.txt");

clarg::argInt plevel     ("-p", "plevel", 1);
clarg::argInt num_it ("-num_it", "number limit of iterations to the CG solver", 500);
clarg::argInt nt_sm("-nt_sm", "Dohr (l1): number of threads to process the submeshes", 0);
clarg::argInt nt_d("-nt_d", "Dohr (l1): number of threads to decompose each submesh", 0);
clarg::argInt nt_a("-nt_a", "Pair (l2): number of threads to assemble each submesh (multiply by nt_sm)", 0);
clarg::argBool dohr_tbb("-dohr_tbb", "assemble TPZDohrStructMatrix (level 1) using TBB", false);
clarg::argBool pair_tbb("-pair_tbb", "assemble TPZPairStructMatrix (level 2) using TBB", false);
clarg::argInt nt_multiply("-nt_m", "number of threads to multiply", 0);
clarg::argInt nsub("-nsub", "number of substructs", 4);

clarg::argInt dim_arg("-dim", "dim???", 3);
clarg::argInt maxlevel("-maxlevel", "maxlevel???", 5);
clarg::argInt sublevel("-sublevel", "sublevel???", 3);

clarg::argInt verb_level("-v", "verbosity level", 0);

clarg::argBool bc("-bc", "binary checkpoints", false);

clarg::argBool h("-h", "help message", false);

clarg::argInt refin("-ref", "refine mesh", 1);

/* Run statistics. */
RunStatsTable total_rst   ("-tot_rdt", "Whole program (total) statistics raw data table");
RunStatsTable create_rst  ("-cre_rdt", "Create statistics raw data table (step 2)");
RunStatsTable assemble_rst("-ass_rdt", "Assemble statistics raw data table (step 3)");
RunStatsTable precond_rst ("-pre_rdt", "Precond statistics raw data table (step 3)");
RunStatsTable solve_rst   ("-sol_rdt", "Solver statistics raw data table (step 4)");

#ifdef USING_LIKWID
#define PERF_START(obj)                                                         \
likwid_markerStartRegion(#obj);                 \
obj.start()
#define PERF_STOP(obj)                                                          \
obj.stop();                                                                     \
likwid_markerStopRegion(#obj)
#else
#define PERF_START(obj) obj.start()
#define PERF_STOP(obj) obj.stop()
#endif


class FileStreamWrapper
{
public:
    FileStreamWrapper() {}
    ~FileStreamWrapper() {}
    
    void OpenWrite(const std::string& fn)
    {
        if (bc.was_set())
            bfs.OpenWrite(fn);
        else
            fs.OpenWrite(fn);
    }
    
    void OpenRead(const std::string& fn)
    {
        if (bc.was_set())
            bfs.OpenRead(fn);
        else
            fs.OpenRead(fn);
    }
    
    operator TPZStream&()
    {
        if (bc.was_set())
            return bfs;
        else
            return fs;
    }
    
protected:
    TPZFileStream  fs;
    TPZBFileStream bfs;
};

#ifdef USING_LIKWID
#include<likwid.h>

struct likwid_manager_t {
    likwid_manager_t() {
        std::cout << "Calling likwid_markerInit()" << std::endl;
        likwid_markerInit();
    }
    ~likwid_manager_t() {
        likwid_markerClose();
        std::cout << "Calling likwid_markerClose()" << std::endl;
    }
};

#endif


int main(int argc, char *argv[])
{
    
    TPZfTime serial;
    
#ifdef USING_LIKWID
    likwid_manager_t likwid_manager;
#endif
    
#ifdef LOG4CXX
    InitializePZLOG("log4cxx.cfg");
#endif
#ifdef USING_TBB
    task_scheduler_init init;
#endif
    
    int main_ret = EXIT_SUCCESS;
    
    /* Parse the arguments */
    if (clarg::parse_arguments(argc, argv)) {
        cerr << "Error when parsing the arguments!" << endl;
        return 1;
    }
    
    if (h.get_value() == true) {
        help(argv[0]);
        return 1;
    }
    
    /* Verbose macro. */
    unsigned verbose = verb_level.get_value();
#   define VERBOSE(level,...) if (level <= verbose) cout << __VA_ARGS__
    
    if (verbose >= 1) {
        std::cout << "- Arguments -----------------------" << std::endl;
        clarg::values(std::cout, false);
        std::cout << "-----------------------------------" << std::endl;
    }
    
    if (!mp.was_set() && !mc.was_set() && !cf1.was_set() &&
        !cf2.was_set() && !cf3.was_set())
    {
        cerr << "A \"starting_point\" must be provided!" << endl;
        help(argv[0]);
        return 1;
    }
    
    PERF_START(total_rst);
    
    if (pair_tbb.was_set())
        TPZPairStructMatrix::gNumThreads = -1;
    else
        TPZPairStructMatrix::gNumThreads = nt_a.get_value();
    
    TPZGeoMesh  *gmesh = 0;
    TPZAutoPointer<TPZCompMesh> cmeshauto = 0;
    TPZDohrStructMatrix* dohrstruct = 0;
    TPZFMatrix<STATE> *rhs = NULL;
    TPZMatrix<STATE> *matptr = 0;
    int dim = dim_arg.get_value();
    TPZCompEl::SetgOrder(plevel.get_value());
    
    bool running = false;
    
    /* Start from malha_cubo or malha_predio? */
    if (mp.was_set() || mc.was_set())
    {
        if (mp.was_set()) // Predio Elastisco
        {
            if (running) {
                cerr << "ERROR: you must select only one of the start modes: "
                << "mp, mc, cf1, cf2 or cf3" << endl;
                exit(1);
            }
            else
                running = true;
            
            gmesh = MalhaPredio();
            cmeshauto = new TPZCompMesh(gmesh);
            cmeshauto->SetDimModel(dim);
            InsertElasticity(cmeshauto);
            cmeshauto->AutoBuild();
        }
        if (mc.was_set()) // Cubo Viscoso
        {
            if (running) {
                cerr << "ERROR: you must select only one of the start modes: "
                << "mp, mc, cf1, cf2 or cf3" << endl;
                exit(1);
            }
            else  running = true;
            
            VERBOSE(1, "Reading MalhaCubo from file: " << mc.get_value() << endl);
            gmesh = MalhaCubo();
            cmeshauto = new TPZCompMesh(gmesh);
            cmeshauto->SetDimModel(dim);
            cmeshauto->SetDefaultOrder(plevel.get_value());
            //cmeshauto->SetAllCreateFunctionsContinuousWithMem();
            //cmeshauto->SetAllCreateFunctionsContinuous();
            InsertViscoElasticityCubo(cmeshauto);
            cmeshauto->AutoBuild();
        }
        
        VERBOSE(1, "Number of equations " << cmeshauto->NEquations() << endl);
        
        dohrstruct = new TPZDohrStructMatrix(cmeshauto);
        dohrstruct->IdentifyExternalConnectIndexes();
        
        VERBOSE(1, "Substructuring the mesh" << endl);
        
        dohrstruct->SubStructure(nsub.get_value());
        
#ifdef LOG4CXX
        {
            std::stringstream str;
            cmeshauto->Print(str);
            LOGPZ_DEBUG(logger,str.str());
        }
#endif
        
        /* Dump checkpoint 1 */
        if (dc1.was_set() && running)
        {
            VERBOSE(1, "Dumping checkpoint 1 into: " << dc1.get_value() << endl);
            FileStreamWrapper CheckPoint1;
            CheckPoint1.OpenWrite(dc1.get_value());
            cmeshauto->Reference()->Write(CheckPoint1, 0);
            cmeshauto->Write(CheckPoint1, 0);
            dohrstruct->Write(CheckPoint1);
        }
        /* Gen/Check checkpoint 1 MD5 signature? */
        if ((gen_sig_ckpt1.was_set() || chk_sig_ckpt1.was_set()) && running)
        {
            TPZMD5Stream sig;
            cmeshauto->Reference()->Write(sig, 0);
            cmeshauto->Write(sig, 0);
            dohrstruct->Write(sig);
            if (chk_sig_ckpt1.was_set()) {
                int ret;
                if ((ret=sig.CheckMD5(chk_sig_ckpt1.get_value()))) {
                    cerr << "ERROR: MD5 Signature for checkpoint 1 does not match. (ret = " << ret << ")" << endl;
                    return 1;
                }
            }
            if (gen_sig_ckpt1.was_set()) {
                int ret;
                if ((ret = sig.WriteMD5(gen_sig_ckpt1.get_value()))) {
                    cerr << "ERROR when writing ckpt 1 MD5 Signature to file (ret = " << ret << "): "
                    << gen_sig_ckpt1.get_value() << endl;
                    return 1;
                }
            }
        }
    }
    
    if(st1.was_set()) running = false;
    
    // Start from Checkpoint 1
    if (cf1.was_set())
    {
        if (running) {
            cerr << "ERROR: you must select only one of the start modes: mp, mc, cf1, cf2 or cf3" << endl;
            exit(1);
        }
        else
            running = true;
        
        gmesh = new TPZGeoMesh;
        cmeshauto = new TPZCompMesh(gmesh);
        dohrstruct = new TPZDohrStructMatrix(cmeshauto);
        /* Read the checkpoint. */
        {
            FileStreamWrapper CheckPoint1;
            CheckPoint1.OpenRead(cf1.get_value());
            gmesh->Read(CheckPoint1, 0);
            cmeshauto->Read(CheckPoint1, gmesh);
            dohrstruct->Read(CheckPoint1);
        }
        
        dim = cmeshauto->Dimension();
        VERBOSE(1, "Reading dim from file. new dim = " << dim << ", old dim = " << dim_arg.get_value() << endl);
        
    }
    
    /* Work between checkpoint 1 and checkpoint 2 */
    if (running) {
        
        PERF_START(create_rst);
        matptr = dohrstruct->Create();
        PERF_STOP(create_rst);
        
        if (dc2.was_set())
        {
            VERBOSE(1, "Dumping checkpoint 2 into: " << dc2.get_value() << endl);
            FileStreamWrapper CheckPoint2;
            CheckPoint2.OpenWrite(dc2.get_value());
            SAVEABLE_STR_NOTE(CheckPoint2,"cmeshauto->Reference()->Write()");
            cmeshauto->Reference()->Write(CheckPoint2, 0);
            SAVEABLE_STR_NOTE(CheckPoint2,"cmeshauto->Write()");
            cmeshauto->Write(CheckPoint2, 0);
            SAVEABLE_STR_NOTE(CheckPoint2,"matptr->Write()");
            matptr->Write(CheckPoint2, 1);
            SAVEABLE_STR_NOTE(CheckPoint2,"dohrstruct->Write()");
            dohrstruct->Write(CheckPoint2);
        }
        
        /* Gen/Check checkpoint 2 MD5 signature? */
        if (gen_sig_ckpt2.was_set() || chk_sig_ckpt2.was_set())
        {
            TPZMD5Stream sig;
            cmeshauto->Reference()->Write(sig, 0);
            cmeshauto->Write(sig, 0);
            matptr->Write(sig, 1);
            dohrstruct->Write(sig);
            
            if (chk_sig_ckpt2.was_set()) {
                if (sig.CheckMD5(chk_sig_ckpt2.get_value())) {
                    cerr << "ERROR: MD5 Signature for checkpoint 2 does not match." << endl;
                    return 1;
                }
            }
            if (gen_sig_ckpt2.was_set()) {
                if (sig.WriteMD5(gen_sig_ckpt2.get_value())) {
                    cerr << "ERROR when writing ckpt 2 MD5 Signature to file: "
                    << gen_sig_ckpt2.get_value() << endl;
                    return 1;
                }
            }
        }
    }
    
    if(st2.was_set()) running = false;
    
    // Start from Checkpoint 2
    if (cf2.was_set())
    {
        if (running) {
            cerr << "ERROR: you must select only one of the start modes: mp, mc, cf1, cf2 or cf3" << endl;
            exit(1);
        }
        else
            running = true;
        
        FileStreamWrapper CheckPoint2;
        CheckPoint2.OpenRead(cf2.get_value());
        gmesh = new TPZGeoMesh;
        SAVEABLE_SKIP_NOTE(CheckPoint2);
        gmesh->Read(CheckPoint2,0);
        cmeshauto = new TPZCompMesh(gmesh);
        SAVEABLE_SKIP_NOTE(CheckPoint2);
        cmeshauto->Read(CheckPoint2, &gmesh);
        SAVEABLE_SKIP_NOTE(CheckPoint2);
        matptr = dynamic_cast<TPZMatrix<STATE> *>(TPZSavable::Restore(CheckPoint2, 0));
        dohrstruct = new TPZDohrStructMatrix(cmeshauto);
        SAVEABLE_SKIP_NOTE(CheckPoint2);
        dohrstruct->Read(CheckPoint2);
    }
    
    cout << "Serial Section: " << serial.ReturnTimeString() << endl;
    
    TPZAutoPointer<TPZMatrix<STATE> > precond = NULL;
    /* Work between checkpoint 2 and checkpoint 3 */
    if (running)
    {
        if (nt_multiply.was_set() && nt_multiply.get_value() != 1)
        {
            dohrstruct->SetNumThreads(nt_multiply.get_value());
        }
        
        PERF_START(assemble_rst);
        TPZAutoPointer<TPZGuiInterface> gui;
        rhs = new TPZFMatrix<STATE>(cmeshauto->NEquations(),1,0.);
        VERBOSE(1,"dohrstruct->Assemble()" << endl);
        if (dohr_tbb.was_set())
            dohrstruct->AssembleTBB(*matptr,*rhs, gui);
        else
            dohrstruct->Assemble(*matptr,*rhs, gui, nt_sm.get_value(), nt_d.get_value());
        PERF_STOP(assemble_rst);
        
        PERF_START(precond_rst);
        precond = dohrstruct->Preconditioner();
        PERF_STOP(precond_rst);
        
        if (dc3.was_set())
        {
            VERBOSE(1, "Dumping checkpoint 3 into: " << dc3.get_value() << endl);
            FileStreamWrapper CheckPoint3;
            CheckPoint3.OpenWrite(dc3.get_value());
            cmeshauto->Reference()->Write(CheckPoint3, 0);
            cmeshauto->Write(CheckPoint3, 0);
            matptr->Write(CheckPoint3, 1);
            precond->Write(CheckPoint3, 1);
            rhs->Write(CheckPoint3, 0);
        }
        
        /* Gen/Check checkpoint 3 MD5 signature? */
        if (gen_sig_ckpt3.was_set() || chk_sig_ckpt3.was_set())
        {
            TPZMD5Stream sig;
            cmeshauto->Reference()->Write(sig, 0);
            cmeshauto->Write(sig, 0);
            matptr->Write(sig, 1);
            precond->Write(sig, 1);
            rhs->Write(sig, 0);
            int ret;
            if (chk_sig_ckpt3.was_set()) {
                if ((ret=sig.CheckMD5(chk_sig_ckpt3.get_value()))) {
                    cerr << "ERROR(ret=" << ret << ") : MD5 Signature for checkpoint 3 does not match." << endl;
                    return 1;
                }
            }
            if (gen_sig_ckpt3.was_set()) {
                if ((ret=sig.WriteMD5(gen_sig_ckpt3.get_value()))) {
                    cerr << "ERROR (ret=" << ret << ") when writing ckpt 3 MD5 Signature to file: "
                    << gen_sig_ckpt3.get_value() << endl;
                    return 1;
                }
            }
        }
    }
    
    if(st3.was_set()) running = false;
    
    // Start from Checkpoint 3
    if (cf3.was_set())
    {
        if (running) {
            cerr << "ERROR: you must select only one of the start modes: mp, mc, cf1, cf2 or cf3" << endl;
            exit(1);
        }
        else
            running = true;
        gmesh = new TPZGeoMesh;
        cmeshauto = new TPZCompMesh(gmesh);
        dohrstruct = new TPZDohrStructMatrix(cmeshauto);
        
        dim = cmeshauto->Dimension();
        VERBOSE(1, "Reading dim from file. new dim = " << dim << ", old dim = " << dim_arg.get_value() << endl);
        
        FileStreamWrapper CheckPoint3;
        CheckPoint3.OpenRead(cf3.get_value());
        gmesh->Read(CheckPoint3, 0);
        cmeshauto->Read(CheckPoint3, gmesh);
        matptr = dynamic_cast<TPZMatrix<STATE> *>(TPZSavable::Restore(CheckPoint3, 0));
        precond = dynamic_cast<TPZMatrix<STATE> *>(TPZSavable::Restore(CheckPoint3, matptr));
        rhs = new TPZFMatrix<STATE>(cmeshauto->NEquations(),1,0.);
        rhs->Read(CheckPoint3, 0);
    }
    
    if (running) {
        
        /* Work after checkpoint 3 */
        TPZAutoPointer<TPZMatrix<STATE> > dohr = matptr;
        
        int neq = dohr->Rows();
        TPZFMatrix<STATE> diag(neq,1,0.), produto(neq,1);
        
        VERBOSE(1, "Number of equations " << neq << endl);
        
        TPZStepSolver<STATE> pre(precond);
        pre.SetMultiply();
        TPZStepSolver<STATE> cg(dohr);
        
        /* Configure the CG solver to iterate:
         - until it converges (residual <= 1.e-8), or
         - until it reaches 500 itearations.
         */
        cg.SetCG(num_it.get_value(),pre,1.e-8,0);
        
        PERF_START(solve_rst);
        cg.Solve(*rhs,diag);
        PERF_STOP(solve_rst);
        
        
        /* checking if the solver converged */
        if (cg.GetTolerance() > 1.e-8)
        {
            cerr << "ERROR: solver do not converged with the limit of iterations."  << endl;
            exit(1);
        }
        
        TPZDohrMatrix<STATE,TPZDohrSubstructCondense<STATE> > *dohrptr =
        dynamic_cast<TPZDohrMatrix<STATE,TPZDohrSubstructCondense<STATE> > *> (dohr.operator->());
        
        if (!dohrptr) {
            DebugStop();
        }
        
        dohrptr->AddInternalSolution(diag);
        
        
        /* Gen/Check checkpoint 4 MD5 signature? */
        if (gen_sig_ckpt4.was_set() || chk_sig_ckpt4.was_set())
        {
            TPZMD5Stream sig;
            
            dohrptr->Write(sig, 0);
            diag.Write(sig, 0);
            cg.Write(sig, 0);
            //cmeshauto->Reference()->Write(sig, 0);
            //cmeshauto->Write(sig, 0);
            //matptr->Write(sig, 1);
            //precond->Write(sig, 1);
            //rhs->Write(sig, 0);
            
            int ret;
            if (chk_sig_ckpt4.was_set()) {
                if ((ret=sig.CheckMD5(chk_sig_ckpt4.get_value()))) {
                    cerr << "ERROR(ret=" << ret << ") : MD5 Signature for checkpoint 4 does not match." << endl;
                    main_ret = 1;
                }
            }
            if (gen_sig_ckpt4.was_set()) {
                if ((ret=sig.WriteMD5(gen_sig_ckpt4.get_value()))) {
                    cerr << "ERROR (ret=" << ret << ") when writting ckpt 4 MD5 Signature to file: "
                    << gen_sig_ckpt4.get_value() << endl;
                    main_ret = 2;
                }
            }
        }
        
        typedef std::list<TPZAutoPointer<TPZDohrSubstructCondense<STATE> > > subtype;
        const subtype &sublist = dohrptr->SubStructures();
        subtype::const_iterator it = sublist.begin();
        int subcount=0;
        while (it != sublist.end()) {
            
            TPZFMatrix<STATE> subext,subu;
            dohrptr->fAssembly->Extract(subcount,diag,subext);
            (*it)->UGlobal(subext,subu);
            TPZCompMesh *submesh = SubMesh(cmeshauto, subcount);
            submesh->LoadSolution(subu);
            
            //ViscoElastico
            //Atualizando memoria do material
            std::map<int ,TPZMaterial * > materialmap(submesh->MaterialVec());
            std::map<int ,TPZMaterial * >::iterator itmat;
            for (itmat = materialmap.begin(); itmat != materialmap.end() ; itmat++)
            {
                TPZMaterial * mat = itmat->second;
                TPZViscoelastic *vmat = dynamic_cast< TPZViscoelastic *> (mat);
                if(vmat)
                {
                    vmat->SetUpdateMem(true);
                }
            }
            subcount++;
            it++;
        }
        
#ifdef LOG4CXX
        {
            std::stringstream sout;
            diag.Print("Resultado do processo iterativo",sout);
            LOGPZ_INFO(loggerconverge,sout.str())
        }
#endif
        
        TPZMaterial * mat = cmeshauto->FindMaterial(1);
        int nstate = mat->NStateVariables();
        int nscal = 0, nvec = 0;
        if(nstate ==1)
        {
            nscal = 1;
        }
        else
        {
            nvec = 1;
        }
        TPZManVector<std::string> scalnames(nscal),vecnames(nvec);
        if(nscal == 1)
        {
            scalnames[0]="state";
        }
        else
        {
            vecnames[0] = "state";
        }
        std::string postprocessname("dohrmann_visco.vtk");
        TPZVTKGraphMesh vtkmesh(cmeshauto.operator->(),dim,mat,scalnames,vecnames);
        vtkmesh.SetFileName(postprocessname);
        vtkmesh.SetResolution(1);
        int numcases = 1;
        
        
        // Iteracoes de tempo
        int istep = 0;
        vtkmesh.DrawMesh(numcases);
        vtkmesh.DrawSolution(istep, 1.);
    }
    
    PERF_STOP(total_rst);
    
    
    if (gmesh != NULL) delete gmesh;
    
    return main_ret;
}

void InsertElasticity(TPZAutoPointer<TPZCompMesh> mesh)
{
                mesh->SetDimModel(3);
                int nummat = 1;
                STATE E = 1.e6;
                STATE poisson = 0.3;
                TPZManVector<STATE> force(3,0.);
                force[1] = 20.;
                TPZElasticity3D *elast = new TPZElasticity3D(nummat,E,poisson,force);
                TPZMaterial * elastauto(elast);
                TPZFMatrix<STATE> val1(3,3,0.),val2(3,1,0.);
                TPZBndCond *bc = elast->CreateBC(elastauto, -1, 0, val1, val2);
                TPZMaterial * bcauto(bc);
                mesh->InsertMaterialObject(elastauto);
                mesh->InsertMaterialObject(bcauto);
}

void InsertViscoElasticity(TPZAutoPointer<TPZCompMesh> mesh)
{
                mesh->SetDimModel(3);
    int nummat = 1;
    TPZManVector<STATE> force(3,0.);
    force[1] = 20.;
    STATE ElaE = 1000., poissonE = 0.2, ElaV = 100., poissonV = 0.1;
    
    STATE lambdaV = 0, muV = 0, alpha = 0, deltaT = 0;
    lambdaV = 11.3636;
    muV = 45.4545;
    alpha = 1.;
    deltaT = 0.01;
    
    TPZViscoelastic *viscoelast = new TPZViscoelastic(nummat);
    viscoelast->SetMaterialDataHooke(ElaE, poissonE, ElaV, poissonV, alpha, deltaT, force);
    
    TPZMaterial * viscoelastauto(viscoelast);
    TPZFMatrix<STATE> val1(3,3,0.),val2(3,1,0.);
    TPZBndCond *bc = viscoelast->CreateBC(viscoelastauto, -1, 0, val1, val2);
    TPZFNMatrix<6,STATE> qsi(6,1,0.);
    viscoelast->SetDefaultMem(qsi); //elast
    viscoelast->PushMemItem(); //elast
    TPZMaterial * bcauto(bc);
    mesh->InsertMaterialObject(viscoelastauto);
    mesh->InsertMaterialObject(bcauto);
}

void InsertViscoElasticityCubo(TPZAutoPointer<TPZCompMesh> mesh)
{
    mesh->SetDimModel(3);
    int nummat = 1, neumann = 1, mixed = 2;
    //      int dirichlet = 0;
    int dir1 = -1, dir2 = -2, dir3 = -3, neumann1 = -4., neumann2 = -5;
    TPZManVector<STATE> force(3,0.);
    //force[1] = 0.;
    REAL Ela = 1000, poisson = 0.;
    REAL lambdaV = 0, muV = 0, alphaT = 0;
    lambdaV = 11.3636;
    muV = 45.4545;
    alphaT = 0.01;
    
    
    //TPZViscoelastic *viscoelast = new TPZViscoelastic(nummat, Ela, poisson, lambdaV, muV, alphaT, force);
    TPZElasticity3D *viscoelast = new TPZElasticity3D(nummat, Ela, poisson, force);
    
    TPZFNMatrix<6> qsi(6,1,0.);
    //viscoelast->SetDefaultMem(qsi); //elast
    //int index = viscoelast->PushMemItem(); //elast
    TPZMaterial * viscoelastauto(viscoelast);
    mesh->InsertMaterialObject(viscoelastauto);
    
    // Neumann em x = 1;
    TPZFMatrix<STATE> val1(3,3,0.),val2(3,1,0.);
    val2(0,0) = 1.;
    TPZBndCond *bc4 = viscoelast->CreateBC(viscoelastauto, neumann1, neumann, val1, val2);
    TPZMaterial * bcauto4(bc4);
    mesh->InsertMaterialObject(bcauto4);
    
    // Neumann em x = -1;
    val2(0,0) = -1.;
    TPZBndCond *bc5 = viscoelast->CreateBC(viscoelastauto, neumann2, neumann, val1, val2);
    TPZMaterial * bcauto5(bc5);
    mesh->InsertMaterialObject(bcauto5);
    
    val2.Zero();
    // Dirichlet em -1 -1 -1 xyz;
    val1(0,0) = 1e4;
    val1(1,1) = 1e4;
    val1(2,2) = 1e4;
    TPZBndCond *bc1 = viscoelast->CreateBC(viscoelastauto, dir1, mixed, val1, val2);
    TPZMaterial * bcauto1(bc1);
    mesh->InsertMaterialObject(bcauto1);
    
    // Dirichlet em 1 -1 -1 yz;
    val1(0,0) = 0.;
    val1(1,1) = 1e4;
    val1(2,2) = 1e4;
    TPZBndCond *bc2 = viscoelast->CreateBC(viscoelastauto, dir2, mixed, val1, val2);
    TPZMaterial * bcauto2(bc2);
    mesh->InsertMaterialObject(bcauto2);
    
    // Dirichlet em 1 1 -1 z;
    val1(0,0) = 0.;
    val1(1,1) = 0.;
    val1(2,2) = 1e4;
    TPZBndCond *bc3 = viscoelast->CreateBC(viscoelastauto, dir3, mixed, val1, val2);
    TPZMaterial * bcauto3(bc3);
    mesh->InsertMaterialObject(bcauto3);
}

TPZGeoMesh *MalhaPredio()
{
                //int nBCs = 1;
                int numnodes=-1;
                int numelements=-1;
                
                string FileName = mp.get_value();
                {
                                bool countnodes = false;
                                bool countelements = false;
                                
                                ifstream read (FileName.c_str());
                                if (!read.is_open()) {
            cerr << "Could not open file: " << FileName << endl;
            exit(1);
                                }
        
                                while(read)
                                {
                                                char buf[1024];
                                                read.getline(buf, 1024);
                                                std::string str(buf);
                                                if(str == "Coordinates") countnodes = true;
                                                if(str == "end coordinates") countnodes = false;
                                                if(countnodes) numnodes++;
                                                
                                                if(str == "Elements") countelements = true;
                                                if(str == "end elements") countelements = false;
                                                if(countelements) numelements++;
                                }
                }
                
                TPZGeoMesh * gMesh = new TPZGeoMesh;
                
                gMesh -> NodeVec().Resize(numnodes);
                
                TPZVec <int64_t> TopolTetra(4);
                
                const int Qnodes = numnodes;
                TPZVec <TPZGeoNode> Node(Qnodes);
                
                //setting nodes coords
                int64_t nodeId = 0, elementId = 0, matElId = 1;
                
                ifstream read;
                read.open(FileName.c_str());
                
                double nodecoordX , nodecoordY , nodecoordZ ;
                
                char buf[1024];
                read.getline(buf, 1024);
                read.getline(buf, 1024);
                std::string str(buf);
                int in;
                for(in=0; in<numnodes; in++)
                {
                                read >> nodeId;
                                read >> nodecoordX;
                                read >> nodecoordY;
                                read >> nodecoordZ;
                                Node[nodeId-1].SetNodeId(nodeId);
                                Node[nodeId-1].SetCoord(0,nodecoordX);
                                Node[nodeId-1].SetCoord(1,nodecoordY);
                                Node[nodeId-1].SetCoord(2,nodecoordZ);
                                gMesh->NodeVec()[nodeId-1] = Node[nodeId-1];
                }
                
                {
                                read.close();
                                read.open(FileName.c_str());
        
                                int l , m = numnodes+5;
                                for(l=0; l<m; l++)
                                {
                                                read.getline(buf, 1024);
                                }
                                
                                int el;
                                int matBCid = -1;
                                //std::set<int> ncoordz; //jeitoCaju
                                for(el=0; el<numelements; el++)
                                {
                                                read >> elementId;
                                                read >> TopolTetra[0]; //node 1
                                                read >> TopolTetra[1]; //node 2
                                                read >> TopolTetra[2]; //node 3
                                                read >> TopolTetra[3]; //node 4
                                                
                                                // O GID comeca com 1 na contagem dos nodes, e nao zero como no PZ, assim o node 1 na verdade é o node 0
                                                TopolTetra[0]--;
                                                TopolTetra[1]--;
                                                TopolTetra[2]--;
                                                TopolTetra[3]--;
                                                
                                                int64_t index = el;
                                                
                                                new TPZGeoElRefPattern< pzgeom::TPZGeoTetrahedra> (index, TopolTetra, matElId, *gMesh);
                                }
                                
                                gMesh->BuildConnectivity();
                                // Colocando as condicoes de contorno
                                
                                for(el=0; el<numelements; el++)
                                {
                                                TPZManVector <TPZGeoNode,4> Nodefinder(4);
                                                TPZManVector <REAL,3> nodecoord(3);
                                                TPZGeoEl *tetra = gMesh->ElementVec()[el];
                                                // na face z = 0
                                                TPZVec<int64_t> ncoordzVec(0); int64_t sizeOfVec = 0;
                                                for (int i = 0; i < 4; i++)
                                                {
                                                                int64_t pos = tetra->NodeIndex(i);
                                                                Nodefinder[i] = gMesh->NodeVec()[pos];
                                                                Nodefinder[i].GetCoordinates(nodecoord);
                                                                if (nodecoord[2] == 0.)
                                                                {
                                                                                sizeOfVec++;
                                                                                ncoordzVec.Resize(sizeOfVec);
                                                                                ncoordzVec[sizeOfVec-1] = pos;
                                                                }
                                                }
                                                if(ncoordzVec.NElements() == 3)
                                                {
                                                                int lado = tetra->WhichSide(ncoordzVec);
                                                                TPZGeoElSide tetraSide(tetra, lado);
                                                                TPZGeoElBC(tetraSide,matBCid);
                                                }
                                }
                }
                
                ofstream arg("malhaPZ.txt");
                gMesh->Print(arg);
                ofstream predio("GeoPredio.vtk");
                TPZVTKGeoMesh::PrintGMeshVTK(gMesh,predio,true);
                
                return gMesh;
}

TPZGeoMesh *MalhaCubo()
{
                int numnodes=-1;
                int numelements=-1;
                
                string FileName = mc.get_value();
                {
                                bool countnodes = false;
                                bool countelements = false;
                                
                                ifstream read (FileName.c_str());
                                if (!read.is_open()) {
            cerr << "Could not open file: " << FileName << endl;
            exit(1);
                                }
                                
                                while(read)
                                {
                                                char buf[1024];
                                                read.getline(buf, 1024);
                                                std::string str(buf);
                                                if(str == "Coordinates") countnodes = true;
                                                if(str == "end coordinates") countnodes = false;
                                                if(countnodes) numnodes++;
                                                
                                                if(str == "Elements") countelements = true;
                                                if(str == "end elements") countelements = false;
                                                if(countelements) numelements++;
                                }
                }
                
                TPZGeoMesh * gMesh = new TPZGeoMesh;
                
                gMesh -> NodeVec().Resize(numnodes);
                
                TPZManVector <int64_t> TopolTetra(4);
                
                const int Qnodes = numnodes;
                TPZVec <TPZGeoNode> Node(Qnodes);
                
                //setting nodes coords
                int nodeId = 0, elementId = 0, matElId = 1;
                
                ifstream read;
                read.open(FileName.c_str());
                
                double nodecoordX , nodecoordY , nodecoordZ ;
                
                char buf[1024];
                read.getline(buf, 1024);
                read.getline(buf, 1024);
                std::string str(buf);
                int in;
                for(in=0; in<numnodes; in++)
                {
                                read >> nodeId;
                                read >> nodecoordX;
                                read >> nodecoordY;
                                read >> nodecoordZ;
                                Node[nodeId-1].SetNodeId(nodeId);
                                Node[nodeId-1].SetCoord(0,nodecoordX);
                                Node[nodeId-1].SetCoord(1,nodecoordY);
                                Node[nodeId-1].SetCoord(2,nodecoordZ);
                                gMesh->NodeVec()[nodeId-1] = Node[nodeId-1];
                }
                
                {
                                read.close();
                                read.open(FileName.c_str());
                                
                                int l , m = numnodes+5;
                                for(l=0; l<m; l++)
                                {
                                                read.getline(buf, 1024);
                                }
                                
                                
                                int el;
                                int neumann1 = -4, neumann2 = -5;
                                //std::set<int> ncoordz; //jeitoCaju
                                for(el=0; el<numelements; el++)
                                {
                                                read >> elementId;
                                                read >> TopolTetra[0]; //node 1
                                                read >> TopolTetra[1]; //node 2
                                                read >> TopolTetra[2]; //node 3
                                                read >> TopolTetra[3]; //node 4
                                                
                                                // O GID comeca com 1 na contagem dos nodes, e nao zero como no PZ, assim o node 1 na verdade é o node 0
                                                TopolTetra[0]--;
                                                TopolTetra[1]--;
                                                TopolTetra[2]--;
                                                TopolTetra[3]--;
                                                
                                                int64_t index = el;
                                                
                                                new TPZGeoElRefPattern< pzgeom::TPZGeoTetrahedra> (index, TopolTetra, matElId, *gMesh);
                                }
                                
                                gMesh->BuildConnectivity();
                                
                                // Colocando as condicoes de contorno
                                for(el=0; el<numelements; el++)
                                {
                                                TPZManVector <TPZGeoNode,4> Nodefinder(4);
                                                TPZManVector <REAL,3> nodecoord(3);
                                                TPZGeoEl *tetra = gMesh->ElementVec()[el];
                                                
                                                // na face x = 1
                                                TPZVec<int64_t> ncoordzVec(0); int64_t sizeOfVec = 0;
                                                for (int i = 0; i < 4; i++)
                                                {
                                                                int64_t pos = tetra->NodeIndex(i);
                                                                Nodefinder[i] = gMesh->NodeVec()[pos];
                                                                Nodefinder[i].GetCoordinates(nodecoord);
                                                                if (nodecoord[0] == 1.)
                                                                {
                                                                                sizeOfVec++;
                                                                                ncoordzVec.Resize(sizeOfVec);
                                                                                ncoordzVec[sizeOfVec-1] = pos;
                                                                }
                                                }
                                                if(ncoordzVec.NElements() == 3)
                                                {
                                                                int lado = tetra->WhichSide(ncoordzVec);
                                                                TPZGeoElSide tetraSide(tetra, lado);
                                                                TPZGeoElBC(tetraSide,neumann1);
                                                }
                                                
                                                // Na face x = -1
                                                ncoordzVec.Resize(0);
                                                sizeOfVec = 0;
                                                for (int i = 0; i < 4; i++)
                                                {
                                                                int pos = tetra->NodeIndex(i);
                                                                Nodefinder[i] = gMesh->NodeVec()[pos];
                                                                
                                                                Nodefinder[i].GetCoordinates(nodecoord);
                                                                if (nodecoord[0] == -1.)
                                                                {
                                                                                sizeOfVec++;
                                                                                ncoordzVec.Resize(sizeOfVec);
                                                                                ncoordzVec[sizeOfVec-1] = pos;
                                                                }
                                                }
                                                if(ncoordzVec.NElements() == 3)
                                                {
                                                                int lado = tetra->WhichSide(ncoordzVec);
                                                                TPZGeoElSide tetraSide(tetra, lado);
                                                                TPZGeoElBC(tetraSide,neumann2);
                                                }
                                                
                                }
                                
                                TPZVec <REAL> xyz(3,-1.), yz(3,-1.), z(3,1.);
                                yz[0] = 1.;
                                z[2] = -1;
                                int bcidxyz = -1, bcidyz = -2, bcidz = -3;
                                SetPointBC(gMesh, xyz, bcidxyz);
                                SetPointBC(gMesh, yz, bcidyz);
                                SetPointBC(gMesh, z, bcidz);
                                
                }
                
                /* refine mesh */
                if (refin.was_set()) {
        
                                int nh = refin.get_value();
        
                                for ( int ref = 0; ref < nh; ref++ ){
                                                TPZVec<TPZGeoEl *> filhos;
                                                int64_t n = gMesh->NElements();
                                                for ( int64_t i = 0; i < n; i++ ){
                                                                TPZGeoEl * gel = gMesh->ElementVec() [i];
                                        
                if(!gel) continue;
                if(gel->Dimension() < 1) continue;
                if(gel->HasSubElement()) continue;
                
                gel->Divide (filhos);
                                                }
                                }
                }
                
                ofstream arg("malhaPZ1BC.txt");
                gMesh->Print(arg);
                
                std::ofstream out("Cube.vtk");
                TPZVTKGeoMesh::PrintGMeshVTK(gMesh, out, true);
                
                return gMesh;
}

void SetPointBC(TPZGeoMesh *gr, TPZVec<REAL> &x, int bc)
{
                // look for an element/corner node whose distance is close to start
                TPZGeoNode *gn1 = gr->FindNode(x);
                int64_t iel;
                int64_t nelem = gr->ElementVec().NElements();
                TPZGeoEl *gel;
                for (iel = 0; iel<nelem; iel++) {
                                gel = gr->ElementVec()[iel];
                                if(!gel) continue;
                                int nc = gel->NCornerNodes();
                                int c;
                                for (c=0; c<nc; c++) {
                                                TPZGeoNode *gn = gel->NodePtr(c);
                                                if (gn == gn1) {
                                                                break;
                                                }
                                }
                                if (c<nc) {
                                                TPZGeoElBC(gel, c, bc);
                                                return;
                                }
                }
}

REAL Height(TPZGeoMesh *gmesh)
{
                TPZAdmChunkVector<TPZGeoNode> &nodevec = gmesh->NodeVec();
                int nnodes = nodevec.NElements();
                int in;
                REAL maxz = 0.;
                for (in=0; in<nnodes; in++) {
                                REAL z = nodevec[in].Coord(2);
                                maxz = (maxz < z) ? z : maxz;
                }
                return maxz;
}

int SubStructure(TPZAutoPointer<TPZCompMesh> cmesh, REAL height)
{
                int64_t nelem = cmesh->NElements();
                TPZManVector<int64_t> subindex(nelem,-1);
                int64_t iel;
                int nsub = 0;
                for (iel=0; iel<nelem; iel++)
                {
                                TPZCompEl *cel = cmesh->ElementVec()[iel];
                                if (!cel) {
                                                continue;
                                }
                                TPZGeoEl *gel = cel->Reference();
                                if (!gel) {
                                                continue;
                                }
                                int nsides = gel->NSides();
                                TPZManVector<REAL> center(gel->Dimension(),0.), xco(3,0.);
                                gel->CenterPoint(nsides-1,center);
                                gel->X(center,xco);
                                REAL z = xco[2];
                                int floor = (int) z/height;
                                nsub = (floor+1) > nsub ? (floor+1) : nsub;
                                subindex[iel] = floor;
                }
                
#ifdef PZDEBUG
                {
                                TPZGeoMesh *gmesh = cmesh->Reference();
                                int64_t nelgeo = gmesh->NElements();
                                TPZVec<int> domaincolor(nelgeo,-999);
                                int64_t cel;
                                int64_t nel = cmesh->NElements();
                                for (cel=0; cel<nel; cel++) {
                                                TPZCompEl *compel = cmesh->ElementVec()[cel];
                                                if(!compel) continue;
                                                TPZGeoEl *gel = compel->Reference();
                                                if (!gel) {
                                                                continue;
                                                }
                                                domaincolor[gel->Index()] = subindex[cel];
                                }
                                ofstream vtkfile("partition.vtk");
                                TPZVTKGeoMesh::PrintGMeshVTK(gmesh, vtkfile, domaincolor);
                }
#endif
                
                int isub;
                TPZManVector<TPZSubCompMesh *> submeshes(nsub,0);
                for (isub=0; isub<nsub; isub++)
                {
                                int64_t index;
                                std::cout << '^'; std::cout.flush();
                                submeshes[isub] = new TPZSubCompMesh(cmesh,index);
                                
                                if (index < subindex.NElements())
                                {
                                                subindex[index] = -1;
                                }
                }
                for (iel=0; iel<nelem; iel++)
                {
                                int domindex = subindex[iel];
                                if (domindex >= 0)
                                {
                                                TPZCompEl *cel = cmesh->ElementVec()[iel];
                                                if (!cel)
                                                {
                                                                continue;
                                                }
                                                submeshes[domindex]->TransferElement(cmesh.operator->(),iel);
                                }
                }
                cmesh->ComputeNodElCon();
                for (isub=0; isub<nsub; isub++)
                {
                                submeshes[isub]->MakeAllInternal();
                                std::cout << '*'; std::cout.flush();
                }
                
                cmesh->ComputeNodElCon();
                cmesh->CleanUpUnconnectedNodes();
                return nsub;
}