output.h

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#ifndef MECHSYS_FEM_OUTPUT_H
#define MECHSYS_FEM_OUTPUT_H

#ifdef HAVE_CONFIG_H
  #include "config.h"
#else
  #ifndef REAL
    #define REAL double
  #endif
#endif

#define OUTPUT_MAX_NODES 20
#define OUTPUT_MAX_ELEMS 10

#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <map>

#include "util/exception.h"
#include "util/numstreams.h"
#include "fem/inputdata.h"
#include "fem/data.h"
#include "fem/solver/intsolverdata.h"
#include "linalg/matrix.h"

using Util::_4;

namespace FEM
{

class Output
{
    friend std::ostream & operator<< (std::ostream & os, FEM::Output const & OU);
public:
    Output() : _n_out_nodes(0), _n_out_elems(0), _notify(NULL), _pt2object(NULL), _pt2wrapper(NULL) {}
    ~Output();
    void Initialize        (FEM::InputData const * idata, FEM::Data const * data);
    void OutputIncrement   (int iStage, int increment);
    void OutputIncrement   (int iStage, int increment, IntSolverData const & ISD);
    void OutputStage       (int iStage);
    void SetNotifyCallBack (void (*pFun)(int, int, IntSolverData const &)) { _notify=pFun; }
    void SetNotifyCallBack (void* pt2Object, void (*pt2Wrapper)(void* pt2Object, int, int) )
    {
         _pt2object = pt2Object;
        _pt2wrapper = pt2Wrapper;
    }
    void OpenFiles();
    void CloseFiles();
private:
    FEM::InputData const * _idata;
    FEM::Data      const * _data;
    int                    _n_out_nodes;
    int                    _n_out_elems;
    std::ofstream          _a_node_files[OUTPUT_MAX_NODES];
    std::ofstream          _a_elem_files[OUTPUT_MAX_ELEMS];
    // VTK
    void VTKCreateFile(String const & Filename);
    void VTKAppendData(String const & Filename, int iStage, int iInc);
    void _write_stage_opendx(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues);
    void _write_stage_vtk   (size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues);
private:
    void (*_notify)(int iStage, int increment, IntSolverData const & ISD);
    void  * _pt2object;
    void (*_pt2wrapper)(void* pt2Object, int iStage, int increment);

}; // class Output

std::ostream & operator<< (std::ostream & os, FEM::Output const & OU)
{
    os << (*OU._data) << std::endl;
    return os;
}




inline void Output::Initialize(FEM::InputData const * idata, FEM::Data const * data) // {{{
{
          _idata = idata;
           _data = data;
    _n_out_nodes = _idata->outNODES.size();
    _n_out_elems = _idata->outELEMS.size();

    // Check number of nodes for output
    if (_n_out_nodes>OUTPUT_MAX_NODES)
        throw new Fatal(_("Output::Initialize: Number of nodes for output too big <max=%d>"),OUTPUT_MAX_NODES);

    // Check range of ID of nodes to output
    for (int i_node=0; i_node<_n_out_nodes; ++i_node)
    {
        int id=_idata->outNODES[i_node]-1;
        if (!(id>=0 && id<static_cast<int>(data->Nodes.size())))
            throw new Fatal(_("Output::Initialize: Number of node < %d > to output is out of range <max=%d>"), id+1, data->Nodes.size());
    }

    // Check number of elements for output
    if (_n_out_elems>OUTPUT_MAX_ELEMS)
        throw new Fatal(_("Output::Initialize: Number of elements for output too big <max=%d>"),OUTPUT_MAX_ELEMS);

    // Check range of ID of elements to output
    for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
    {
        int id=_idata->outELEMS[i_elem]-1;
        if (!(id>=0 && id<static_cast<int>(data->Elements.size())))
            throw new Fatal(_("Output::Initialize: Number of element < %d > to output is out of range <max=%d>"), id+1, data->Elements.size());
    }
    //OpenFiles(); // 

} // }}}

Output::~Output() // {{{
{
    CloseFiles();
} // }}}

void Output::OutputIncrement(int iStage, int increment) // {{{
{
    IntSolverData isd;
    OutputIncrement(iStage, increment, isd);
} // }}}

void Output::OutputIncrement(int iStage, int increment, IntSolverData const & ISD) // {{{
{
    // Write VTK file
    if (increment==-1) VTKCreateFile(_idata->fnOutVTK);
    else
    {
        if (_idata->vtkOutAllIncs) VTKAppendData(_idata->fnOutVTK, iStage, increment);
        else
            if (increment==_idata->nDIV[iStage]-1)
                VTKAppendData(_idata->fnOutVTK, iStage, increment);
    }

    // Loop along nodes for output
    for (int i_node=0; i_node<_n_out_nodes; ++i_node)
    {
        if (_data->Nodes[_idata->outNODES[i_node]-1].Refs()>0)
        {
            if (increment==-1)
                _a_node_files[i_node] << _6()<< "Inc"       << _data->Nodes[ _idata->outNODES[i_node]-1 ].Out(true);
                _a_node_files[i_node] << _6()<< increment+1 << _data->Nodes[ _idata->outNODES[i_node]-1 ].Out();
        }
    }

    // Loop along elements for output
    for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
    {
        if (_data->Elements[_idata->outELEMS[i_elem]-1]->IsActive())
        {
            if (increment==-1)
                _a_elem_files[i_elem] << _6()<< "Inc"       << _data->Elements[ _idata->outELEMS[i_elem]-1 ]->OutCenter(true);
                _a_elem_files[i_elem] << _6()<< increment+1 << _data->Elements[ _idata->outELEMS[i_elem]-1 ]->OutCenter();
        }
    }
    // Notify
    if (_notify   !=NULL) (*_notify)(iStage, increment, ISD);
    if (_pt2object!=NULL) { if (_pt2wrapper!=NULL) (*_pt2wrapper)(_pt2object, iStage, increment); }

} // }}}

void Output::OutputStage(int iStage) // {{{
{
    if (_idata->DoOutSTAGE[iStage]==0) return;

    int n_nodes  = _data->Nodes.size();
    Array<Element *> const & elements = _data->Elements;
    Array<Element *> active_elements;
    int n_elements                   = elements.size();
    int n_active_elems               = 0;
    
    // Filter elements
    // Inactive elements are not considered!
    // Line and Embedded type elements do not be considered too!
    for (int i_elem=0; i_elem<n_elements; ++i_elem)
    {
        String ele_name = elements[i_elem]->Name();
        if (elements[i_elem]->IsActive() == false) continue;
        if (ele_name.substr(0,3)         == "Lin") continue;
        active_elements.push_back(elements[i_elem]);
        n_active_elems++;
    }

    int const MAX_DATA_COMP = 20;
    LinAlg::Matrix<REAL  > values(n_nodes, MAX_DATA_COMP);
    LinAlg::Matrix<size_t> refs  (n_nodes, MAX_DATA_COMP);
    values.SetValues(0.0);
    values.SetValues(0  );

    std::map<String, int>  index_map;
    Array<String>          all_labels;
    int total_comps = 0;   // number of total components for all element types used

    // Loop in elements to get nodal information
    for (int i_elem=0; i_elem<n_active_elems; ++i_elem)
    {
        LinAlg::Matrix<REAL>    elem_values;
        Array<String>           elem_labels;
        active_elements[i_elem]->OutNodes(elem_values, elem_labels);
        int n_labels                  = elem_labels.size();
        Array<Node*> const & connects = active_elements[i_elem]->Connects();
        int              n_elem_nodes = active_elements[i_elem]->nNodes();
        for (int j_label = 0; j_label < n_labels; ++j_label)
        {
            String & current_label = elem_labels[j_label];
            if (index_map.find(elem_labels[j_label]) == index_map.end()) 
            {
                total_comps++;
                all_labels.push_back(elem_labels[j_label]);
                index_map[current_label] = total_comps-1;
            }

            int index = index_map[current_label];
            
            for (int j_node=0; j_node<n_elem_nodes; ++j_node)
            {
                int node_number            = connects[j_node]->Number();
                values(node_number,index) += elem_values(j_node, j_label);
                refs  (node_number,index) ++;
            }
        }
    }
    assert(total_comps<MAX_DATA_COMP);
    
    //Averaging values
    for (int i=0; i<n_nodes; i++)
        for (int j=0; j<total_comps; j++)
        {
            if (refs(i,j) != 0)
                values(i,j) /= refs(i,j);
            else
                values(i,j)  = 0.0; 
        }
    
    _write_stage_vtk   (iStage, all_labels, values);
    _write_stage_opendx(iStage, all_labels, values);

} // }}}

inline void Output::_write_stage_opendx(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues) // {{{
{
    // Open/create file
    std::ofstream ofile;
    std::ostringstream oss;
    oss << _idata->fnOutSTAGES[iStage].GetSTL().c_str() << ".dx";
    ofile.open(oss.str().c_str(), std::ios::out);

    // Arrays of nodes and elements
    Array<Node>     const &     N = _data->Nodes;
    Array<Element*> const & All_E = _data->Elements; //all elements with out filtering
    Array<Element*> E;

    // Number of nodes and elements and data components
    size_t n_nodes     = N.size();
    size_t n_all_elems = All_E.size();
    size_t n_elems     = 0; // only active elements
    size_t n_comps     = Labels.size(); // number of components
    
    // Filter elements
    // Inactive elements are not considered!
    // Line and Embedded type elements do not be considered too!
    for (size_t i=0; i<n_all_elems; ++i)
    {
        String ele_name =All_E[i]->Name();
        if (All_E[i]->IsActive() == false) continue;
        if (ele_name.substr(0,3)  == "Lin") continue;
        E.push_back(All_E[i]);
        n_elems++;
    }
    
    // Write opendx file
    // Writing Nodes
    // opendx can only handle one type of element at the same type.
    // there is no shape for elements with 20 nodes.
                                                     ofile << "object \"nodes\" class array type float rank 1 shape 3 items " << n_nodes << " data follows"       << std::endl;
    for (size_t i=0; i<n_nodes; ++i)               { ofile << _8s()<<N[i].X() << _8s() << N[i].Y() << _8s() << N[i].Z()                                           << std::endl; }
                                                     ofile <<                                                                                                        std::endl;
                                                     ofile << "object \"elements\" class array type int rank 1 shape ";
    if (E[0]->nNodes() != 20)                        ofile << E[0]->nNodes();
    else                                             ofile << 8;
                                                     ofile << " items " << n_elems << " data follows"                                                      << std::endl;
    for (size_t i=0; i<n_elems; i++)               { 
        size_t n_nodes = E[i]->nNodes();             // reorganize and write connectivities {{{
        Array<Node*> const & connects = \
        E[i]->Connects();
        if (n_nodes==8) {
                                                     ofile << _6() << connects[0]->Number(); //numeration in opendx starts at zero.
                                                     ofile << _6() << connects[1]->Number();
                                                     ofile << _6() << connects[3]->Number();
                                                     ofile << _6() << connects[2]->Number();
                                                     ofile << _6() << connects[4]->Number();
                                                     ofile << _6() << connects[5]->Number();
                                                     ofile << _6() << connects[7]->Number();
                                                     ofile << _6() << connects[6]->Number(); }
        else if (n_nodes==20) {
                                                     ofile << _6() << connects[0] ->Number();
                                                     ofile << _6() << connects[2] ->Number();
                                                     ofile << _6() << connects[6] ->Number();
                                                     ofile << _6() << connects[4] ->Number();
                                                     ofile << _6() << connects[12]->Number();
                                                     ofile << _6() << connects[14]->Number();
                                                     ofile << _6() << connects[18]->Number();
                                                     ofile << _6() << connects[16]->Number(); }
        else { 
            for (int j=0; j<E[i]->nNodes(); ++j)   { ofile << E[i]->Connects()[j]->Number() << " ";                                                                             } } 
                                                     ofile << std::endl;                                                                                                        } 
                                                     ofile << "attribute \"element type\" string \"cubes\"" <<                                                       std::endl;
                                                     ofile << "attribute \"ref\" string \"positions\"" <<                                                            std::endl;
                                                     ofile <<                                                                                                        std::endl;
    for (size_t i=0; i<n_comps; i++)               { ofile << "object \"" << Labels[i] << "\" class array type float rank 0 items " << n_nodes << " data follows" << std::endl;
        for (size_t j=0; j<n_nodes; ++j)           { ofile << _8s() << DataValues(j,i)                                                                            << std::endl; } }
                                                     ofile << std::endl;
                                                     ofile << "object \"_field\" class field"                                                                     << std::endl;
                                                     ofile << "component \"positions\" value \"nodes\""                                                           << std::endl;  
                                                     ofile << "component \"connections\" value \"elements\""                                                      << std::endl;
                                                     ofile << "component \"data\" value \"" << Labels[0] << "\""                                                  << std::endl;
    // Close file
    ofile.close();

} // }}}

inline void Output::_write_stage_vtk(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues) // {{{
{ 
    // Open/create file
    std::ofstream ofile;
    std::ostringstream oss;
    oss << _idata->fnOutSTAGES[iStage].GetSTL().c_str() << ".vtk";
    ofile.open(oss.str().c_str(), std::ios::out);

    // Arrays of nodes and elements
    Array<Node>     const & N = _data->Nodes;
    Array<Element*> const & E = _data->Elements;

    // Number of nodes and elements and data components
    size_t n_nodes = N.size();
    size_t n_all_elems = E.size();
    size_t n_elems = 0;
    size_t n_comps = Labels.size();
    
    // Filter elements
    Array<Element * > active_E;
    for (size_t i=0; i <n_all_elems; ++i)
    {
        String ele_name = E[i]->Name();
        if (E[i]->IsActive()     == false) continue; // Inactive elements are not considered!
        if (ele_name.substr(0,3) == "Lin") continue; // Line and Embedded type elements do not be considered too!
        active_E.push_back(E[i]);
        n_elems++;
    }
    
    // Total number of CELLS data = Sum {1 + nNodes}; 1 for the numPts label
    int n_data = 0;
    for (size_t i=0; i<n_elems; ++i) n_data += 1 + active_E[i]->nNodes();

    // Defining variables for displacements
    const String DUX = "Dux";
    const String DUY = "Duy";
    const String DUZ = "Duz";

    // Write VTK file
                                                     ofile << "# vtk DataFile Version 3.0"                              << std::endl;
                                                     ofile << "MechSys/FEM - " << _idata->projNAME                      << std::endl;
                                                     ofile << "ASCII"                                                   << std::endl;
                                                     ofile << "DATASET UNSTRUCTURED_GRID"                               << std::endl;
                                                     ofile <<                                                              std::endl;
                                                     ofile << "POINTS " << n_nodes << " float"                          << std::endl;
    for (size_t i=0; i<n_nodes; ++i)               { ofile << _8s()<<N[i].X() << _8s() << N[i].Y() << _8s() << N[i].Z() << std::endl; }
                                                     ofile <<                                                              std::endl;
                                                     ofile << "CELLS "<< n_elems << " " << n_data                       << std::endl; 
    for (size_t i=0; i<n_elems; ++i)               { ofile << active_E[i]->nNodes() << " ";  
    for (int    j=0; j<active_E[i]->nNodes(); ++j) { ofile << active_E[i]->Connects()[j]->Number() << " ";                            }         
                                                     ofile <<                                                              std::endl; }
                                                     ofile <<                                                              std::endl;
                                                     ofile << "CELL_TYPES " << n_elems                                  << std::endl; 
    for (size_t i=0; i<n_elems; ++i)               { ofile << active_E[i]->VTKCellType()                                << std::endl; }
                                                     ofile <<                                                              std::endl;
                                                     ofile << "POINT_DATA " << n_nodes                                  << std::endl;

                                                   { ofile << "VECTORS "    << "Disp float"                             << std::endl;
        for (size_t j=0; j<n_nodes; ++j)
            if (N[j].HasVar(DUX) && N[j].HasVar(DUY) && N[j].HasVar(DUZ))
                                                   { ofile << _8s() << N[j].DOFVar(DUX).EssentialVal   
                                                           << _8s() << N[j].DOFVar(DUY).EssentialVal
                                                           << _8s() << N[j].DOFVar(DUZ).EssentialVal                    << std::endl; }
            else                                   { ofile << _8s() << 0 << _8s() << 0 << _8s() << 0                    << std::endl;  }
                                                     ofile <<                                                              std::endl; }

    for (size_t i=0; i<n_comps; i++)               { ofile << "SCALARS "    << Labels[i] << " float 1"                  << std::endl;
                                                     ofile << "LOOKUP_TABLE default"                                    << std::endl; 
        for (size_t j=0; j<n_nodes; ++j)           { ofile << _8s() << DataValues(j,i)                                  << std::endl; }
                                                     ofile <<                                                              std::endl; }
    // Close file
    ofile.close();

} // }}}

inline void Output::OpenFiles() // {{{
{
    // Loop along nodes for output and open respective files
    for (int i_node=0; i_node<_n_out_nodes; ++i_node)
    {
        // Open/create file
        _a_node_files[i_node].open(_idata->fnOutNODES[i_node].GetSTL().c_str(), std::ios::out);
        if (_a_node_files[i_node].fail())
            throw new Fatal(_("Output::Initialize: Could not open file < %s >"), _idata->fnOutNODES[i_node].c_str());
    }

    // Loop along elements for output and open respective files
    for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
    {
        // Open/create file
        _a_elem_files[i_elem].open(_idata->fnOutELEMS[i_elem].GetSTL().c_str(), std::ios::out);
        if (_a_elem_files[i_elem].fail())
            throw new Fatal(_("Output::Initialize: Could not open file < %s >"), _idata->fnOutELEMS[i_elem].c_str());
    }
} // }}}

inline void Output::CloseFiles() // {{{
{
    // Loop along nodes for output and close respective files
    for (int i_node=0; i_node<_n_out_nodes; ++i_node)
        _a_node_files[i_node].close();

    // Loop along elements for output and close respective files
    for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
        _a_elem_files[i_elem].close();
} // }}}

inline void Output::VTKCreateFile(String const & Filename) // {{{
{
    // Open file
    std::ofstream of(Filename.GetSTL().c_str(), std::ios::out);

    // Arrays of nodes and elements
    Array<FEM::Node>     const & N = _data->Nodes;
    Array<FEM::Element*> const & E = _data->Elements;

    // Number of nodes and elements
    int n_nodes = N.size();
    int n_elems = E.size();

    // Total number of CELLS data = Sum {1 + nNodes}; 1 for the numPts label
    int n_data = 0;
    for (int i=0; i<n_elems; ++i) n_data += 1 + E[i]->nNodes();

    // Write VTK file
    of << "# vtk DataFile Version 3.0"                            << "\n";
    of << "MechSys/FEM - "<<_idata->projNAME                      << "\n";
    of << "ASCII"                                                 << "\n";
    of << "DATASET UNSTRUCTURED_GRID"                             << "\n";
    of <<                                                            "\n";
    of << "POINTS "<<n_nodes<<" float"                            << "\n"; for (int i=0; i<n_nodes; ++i) {
    of << _8s()<<N[i].X() <<_8s()<<N[i].Y() <<_8s()<<N[i].Z()     << "\n"; }
    of <<                                                            "\n";
    of << "CELLS "<<n_elems<<" "<<n_data                          << "\n"; for (int i=0; i<n_elems; ++i) {
    of << E[i]->nNodes()<< " ";  for (int j=0; j<E[i]->nNodes(); ++j)
          { of<<E[i]->Connects()[j]->Number()<< " "; }         of << "\n"; }
    of <<                                                            "\n";
    of << "CELL_TYPES " << n_elems                                << "\n"; for (int i=0; i<n_elems; ++i) {
    of << E[i]->VTKCellType()                                     << "\n"; }
    of <<                                                            "\n";
    of << "POINT_DATA " << n_nodes                                << "\n";
    of << "SCALARS Zcoordinate float 1"                           << "\n";
    of << "LOOKUP_TABLE default"                                  << "\n"; for (int i=0; i<n_nodes; ++i) {
    of << _8s()<<N[i].Z()                                         << "\n"; }

    // Close file
    of.close();

} // }}}

inline void Output::VTKAppendData(String const & Filename, int iStage, int iInc) // {{{
{
    // Open file
    std::ofstream of(Filename.GetSTL().c_str(), std::ios::app);

    // Arrays of nodes and elements
    //Array<FEM::Node>     const & N = _data->Nodes;
    Array<FEM::Element*> const & E = _data->Elements;

    // Number of nodes and elements
    //int n_nodes = N.size();
    int n_elems = E.size();

    // Append to VTK file
    of <<                                                             "\n";
    of << "CELL_DATA " << n_elems                                  << "\n";     if (_idata->vtkOutScalar1) {
    of << "SCALARS "<<iStage<<"-"<<iInc<<"-Scalar1 float 1"        << "\n";
    of << "LOOKUP_TABLE default"                                   << "\n"; for (int i=0; i<n_elems; ++i) {
    of << _8s()<<E[i]->OutScalar1()                                << "\n"; } } if (_idata->vtkOutScalar2) {
    of << "SCALARS "<<iStage<<"-"<<iInc<<"-Scalar2 float 1"        << "\n";
    of << "LOOKUP_TABLE default"                                   << "\n"; for (int i=0; i<n_elems; ++i) {
    of << _8s()<<E[i]->OutScalar2()                                << "\n"; } } if (_idata->vtkOutTensor1) {
    of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor1 float"          << "\n"; for (int i=0; i<n_elems; ++i) {
    String ten; E[i]->OutTensor1(ten);                   of << ten << "\n"; } } if (_idata->vtkOutTensor2) {
    of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor2 float"          << "\n"; for (int i=0; i<n_elems; ++i) {
    String ten; E[i]->OutTensor2(ten);                   of << ten << "\n"; } } if (_idata->vtkOutTensor3) {
    of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor3 float"          << "\n"; for (int i=0; i<n_elems; ++i) {
    String ten; E[i]->OutTensor3(ten);                   of << ten << "\n"; } }
    of <<                                                             "\n";
    //of << "POINT_DATA " << n_nodes                                 << "\n";     if (_idata->vtkOutEssen) {
    //of << "VECTORS "<<iStage<<"-"<<iInc<<"-EssentialVectors float" << "\n"; for (int i=0; i<n_nodes; ++i) {
    //String vec; N[i].OutEssentialVector(vec);            of << vec << "\n"; } }  if (_idata->vtkOutNatur) {
    //of << "VECTORS "<<iStage<<"-"<<iInc<<"-NaturalVectors float"   << "\n"; for (int i=0; i<n_nodes; ++i) {
    //String vec; N[i].OutNaturalVector(vec);              of << vec << "\n"; } }

    // Close file
    of.close();

} // }}}

}; // namespace FEM

#endif // MECHSYS_FEM_OUTPUT_H

// vim:fdm=marker

---