mesh.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           mesh.cpp
// author(s):      Ladislav Svoboda
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program 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 Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

// toto je kvuli WIN32_LEAN_AND_MEAN
#include "types.h"

//included libraries
#include "mesh.h"

#include "problem.h"
#include "mnode.h"
#include "melement.h"
#include "vtk.h"
#include "esuf.h"
#include "matrixOperations.h"

#include "librw.h"
#include "gelib.h"

#include <string.h>


namespace mumech {

Mesh :: Mesh (long i, const Problems *p)
{
  // ASSIGNED PROBLEMS
  id = i;
  npa = 5;
  np = 0;
  P = new const Problems*[npa];  memset (P, 0, npa*sizeof(Problem*));
  
  // MESH GEOMETRY
  coll = -1;

  md = 0;
  mtA = MT_VOID;

  volume = -1;

  // MESH COMPONENTS
  nNodes = nElems = 0;
  Nodes = NULL;
  Elems = NULL;

  // MESH RESULTS
  nLC = -77;
  no_dspl_rf = no_strn_rf = no_strs_rf = NULL;
  el_dspl_rf = el_strn_rf = el_strs_rf = NULL;
  
  //displc = strain = stress = NULL;
  
  //
  if (p)
    this->add_problem(p);
  
}
Mesh :: Mesh (long i, const Problems *p, const Mesh *src)
{
  // ASSIGNED PROBLEMS
  id = i;
  npa = 5;
  np = 0;
  P = new const Problems*[npa];  memset (P, 0, npa*sizeof(Problems*));

  // MESH GEOMETRY
  coll = src->coll;

  md = src->md;
  mtA = src->mtA;    if(mtA != MT_GENERAL) _errorr("only general mesh supported in copy constructor");

  volume = -1;
  
  // MESH COMPONENTS
  nNodes = src->nNodes;
  Nodes = new mNode*[nNodes];
  for (long i=0; i<nNodes; i++)
    Nodes[i] = new mNode(i, this, src->Nodes[i]);

  nElems = src->nElems;
  Elems = new mElement*[nElems];
  for (long i=0; i<nElems; i++)
    Elems[i] = new mElement(i, this, src->Elems[i]);
  
  // MESH RESULTS
  nLC = -77;   if (src->nLC != -77)  errol;
  no_dspl_rf = no_strn_rf = no_strs_rf = NULL;   if (src->no_dspl_rf != NULL  ||  src->no_strn_rf != NULL  ||  src->no_strs_rf != NULL)  _errorr("results data not supported in copy constructor");
  el_dspl_rf = el_strn_rf = el_strs_rf = NULL;   if (src->el_dspl_rf != NULL  ||  src->el_strn_rf != NULL  ||  src->el_strs_rf != NULL)  _errorr("results data not supported in copy constructor");
  
  //displc = strain = stress = NULL;
  
  //
  if (p)
    this->add_problem(p);
  
}

Mesh :: ~Mesh ()
{
  memset (P, 0, npa*sizeof(Problem*));
  delete [] P;  
  
  if (Nodes)  for (int i=0; i<nNodes; i++)  delete Nodes[i];
  if (Elems)  for (int i=0; i<nElems; i++)  delete Elems[i];
  
  delete [] Nodes;
  delete [] Elems;
  
  DeleteArray2D (no_dspl_rf, npa);
  DeleteArray2D (no_strn_rf, npa);
  DeleteArray2D (no_strs_rf, npa);
  DeleteArray2D (el_dspl_rf, npa);
  DeleteArray2D (el_strn_rf, npa);
  DeleteArray2D (el_strs_rf, npa);
  
  // if (displc) { for (int i=0; i<nLC; i++)  delete [] displc[i];  delete [] displc; }
  // if (strain) { for (int i=0; i<nLC; i++)  delete [] strain[i];  delete [] strain; }
  // if (stress) { for (int i=0; i<nLC; i++)  delete [] stress[i];  delete [] stress; }
}

int Mesh :: add_problem (const Problems *p)
{
  if (p->is_converted_to_equivalent() == false)  errol;
  
  // number of load cases
  if (this->nLC == -77)
    this->nLC = p->give_nLC();
  else
    if (this->nLC != p->give_nLC())
      errol;
  
  
  //
  P[np++] = p;
  return np-1;
}

//
void Mesh :: set_coll (int val)
{
  if (coll != -1)  errol;
  coll = val;

  if (md != 0)  errol;
  if      (coll == 3)  md = 2;
  else if (coll == 0)  md = 3;
  else errol;
}



bool Mesh :: equal_dimensions (int pid) const
{
  return  this->md == this->P[pid]->give_dimension();
}


bool Mesh :: is_unset (void)
{
  return ! ( coll != -1  ||  mtA != MT_VOID  ||
         volume != -1.0  ||
         nNodes != 0  ||  nElems != 0  ||
         Nodes != NULL  ||  Elems != NULL  ||
         no_dspl_rf != NULL  ||  no_strn_rf != NULL  ||  no_strs_rf != NULL  ||
         el_dspl_rf != NULL  ||  el_strn_rf != NULL  ||  el_strs_rf != NULL
         );
}

void Mesh :: scale (const double *s)
{
  bool twd = this->md == 2;
  
  for (long i=0; i<nNodes; i++) {
    Nodes[i]->coords.x *= s[0];
    Nodes[i]->coords.y *= s[1];
    if (twd) Nodes[i]->coords.z = 0.0;
    else     Nodes[i]->coords.z *= s[2];
  }
}

void Mesh :: local2global (const double *origin, const double *TInv)
{
  bool twd = this->md == 2;
  
  for (long i=0; i<nNodes; i++) {
    // rotating of the coordinate system
    if (twd) Nodes[i]->coords.beRotatedBy2d(TInv);
    else     Nodes[i]->coords.beRotatedBy(TInv);
    
    // swap of coordinate systems
    Nodes[i]->coords.add(origin);
  }
}

void Mesh :: set_and_alloc_nNodes (long nn)
{
  nNodes = nn;
  if (Nodes) errol;
  Nodes = new mNode*[nNodes];
  for (long i=0; i<nNodes; i++)
    Nodes[i] = new mNode(i, this);
}
void Mesh :: set_and_alloc_nElems (long ne)
{
  nElems = ne;
  if (Elems) errol;
  Elems = new mElement*[nElems];
  for (long i=0; i<nElems; i++)
    Elems[i] = new mElement(i, this);
}

void Mesh :: shift_id (long snn, long sne)
{
  for (long i=0; i<nNodes; i++)
    Nodes[i]->id += snn;

  for (long i=0; i<nElems; i++) {
    Elems[i]->id += sne;
    AddScalarToVector(snn, Elems[i]->nodes, Elems[i]->nnodes);
  }
}
  

void Mesh :: generate_regular_mesh (const double *p1, const double *p2, const long *n)
{
  if (this->geom_is_readed()) _errorr("mesh data already given");
  mtA = MT_REGULAR;
  
  ;    if ( n[0] &&  n[1] &&  n[2])  this->set_coll(0);
  else if (!n[0] &&  n[1] &&  n[2])  this->set_coll(1);
  else if ( n[0] && !n[1] &&  n[2])  this->set_coll(2);
  else if ( n[0] &&  n[1] && !n[2])  this->set_coll(3);
  else _errorr("");
  
  long i, j, k, l;
  
  //
  CopyVector (p1, point1, 3);
  CopyVector (p2, point2, 3);
  CopyVector (n,  nsegs,  3);
  if (coll == 1) delta[0] = 0.0;  else  delta[0] = (p2[0] - p1[0]) / n[0];
  if (coll == 2) delta[1] = 0.0;  else  delta[1] = (p2[1] - p1[1]) / n[1];
  if (coll == 3) delta[2] = 0.0;  else  delta[2] = (p2[2] - p1[2]) / n[2];
  
  
  // NODES
  this->set_and_alloc_nNodes ( (n[0] + 1) * (n[1] + 1) * (n[2] + 1) );
  
  //
  l = 0;
  for (k=0; k<n[2]+1; k++)
    for (j=0; j<n[1]+1; j++)
      for (i=0; i<n[0]+1; i++) {
    Nodes[l]->coords[0] = p1[0] + i * delta[0];
    Nodes[l]->coords[1] = p1[1] + j * delta[1];
    Nodes[l]->coords[2] = p1[2] + k * delta[2];
    l++;
      }
  
  
  // ELEMENTS
  this->set_and_alloc_nElems( (n[0] ? n[0] : 1) * (n[1] ? n[1] : 1) * (n[2] ? n[2] : 1) );
  
  // element type
  int et;
  if (coll)  et = 9;   // 2D -  9 - quad       (8  - pixel, ma otoceny cislovani uzlu!!)
  else        et = 12;  // 3D - 12 - hexahedron (11 - voxel, ma otoceny cislovani uzlu!!)
  //
  for (i=0; i<nElems; i++)
    Elems[i]->set_type(et);
  
  // element nodes
  l = 0;
  if (coll == 0) {
    for (i=0; i<n[0]; i++)
      for (j=0; j<n[1]; j++)
    for (k=0; k<n[2]; k++) {
      Elems[l]->nodes[0] = i   +  j   *(n[0]+1) +  k   *(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[1] = i+1 +  j   *(n[0]+1) +  k   *(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[2] = i+1 + (j+1)*(n[0]+1) +  k   *(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[3] = i   + (j+1)*(n[0]+1) +  k   *(n[0]+1)*(n[1]+1) ;
      
      Elems[l]->nodes[4] = i   +  j   *(n[0]+1) + (k+1)*(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[5] = i+1 +  j   *(n[0]+1) + (k+1)*(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[6] = i+1 + (j+1)*(n[0]+1) + (k+1)*(n[0]+1)*(n[1]+1) ;
      Elems[l]->nodes[7] = i   + (j+1)*(n[0]+1) + (k+1)*(n[0]+1)*(n[1]+1) ;
      l++;
    }
  }
  else if (coll == 3) {
    for (i=0; i<n[0]; i++)
      for (j=0; j<n[1]; j++) {
    Elems[l]->nodes[0] = i   +  j   *(n[0]+1) ;
    Elems[l]->nodes[1] = i+1 +  j   *(n[0]+1) ;
    Elems[l]->nodes[2] = i+1 + (j+1)*(n[0]+1) ;
    Elems[l]->nodes[3] = i   + (j+1)*(n[0]+1) ;
    l++;
      }
  }
  else _errorr("");
  
}

int Mesh :: file_type_s2e (const char *s)
{
  if      (_STRCMP(s, "2D"   ) == 0)  return 3;
  else if (_STRCMP(s, "3D"   ) == 0)  return 0;
  else _errorr2 ("Invalid structure of the given VTK file, the key word %s at the second line is not supported", s);
  return -1;
}



enum MeshTypeRead { MTR_VOID, MTR_geom, MTR_fem, MTR_femadd };

void Mesh :: read_geometry_file_vtk (const char *vtkMeshFile, int pid, int lc, bool add)
{
  if (this->geom_is_readed()) {
    if (!add)  _errorr("mesh data already readed");
  }
  else {
    if (add)  _errorr("mesh data not readed");
    else      mtA = MT_GENERAL;
  }
  
  // reading type
  MeshTypeRead mtr = MTR_VOID;
  if      (pid == -1  &&  lc == -1  &&  add == false)
    mtr = MTR_geom;
  else if (pid != -1  &&  lc != -1  &&  add == false)
    mtr = MTR_fem;
  else if (pid != -1  &&  lc != -1  &&  add == true)
    mtr = MTR_femadd;
  else
    errol;
  
  
  //if (PP[0]->give_verbose() > 0)
  //  printf( "Reading of VTK file of points, please wait ...\n" );
  
  // open stream with automatic detection VTK legaci/XML file format
  char *filename = new char[255];
  strcpy(filename, vtkMeshFile);
  Stream *stream = new Stream();
  stream->open(STRM_void, "r", (const char*&)filename); // do not delete filename, it is deleted inside of the function
  
  //
  char *WORD=NULL;
  const XMLElement *xnel = NULL;
  
  bool lgc = stream->isFile();
  
  if (lgc) WORD = new char[255];
  
  
  // HEAD OD FILE
  if (lgc) {
    // 1. line
    FP_skip_line (stream->file());
    
    // 2. line
    // Only strictly 2d or 3d meshes are supported. So I want to know it in advance by keyword in comment line.
    FP_scan_word (stream->file(), WORD);
    FP_skip_line (stream->file());
    
    int cll = this->file_type_s2e (WORD);
    if (!add)
      this->set_coll (cll);
    else
      if (cll != this->coll)  errol;
        
    // 3. line
    FP_scan_expected_word_exit (stream->file(), "ASCII", "Invalid structure of the given VTK file", true);
    FP_skip_line (stream->file());
    // 4. line
    FP_scan_expected_word_exit (stream->file(), "DATASET UNSTRUCTURED_GRID", "Invalid structure of the given VTK file", true);
    FP_skip_line (stream->file());
  }
  else {
    // check head of xml file and read possible ctrl section
    // return XMLElement "Piece"
    const XMLNode *parent = stream->tixnod();
    
    //* declaration - delete if present
    if (parent->FirstChild()->ToDeclaration() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
    
    //* comment - delete if present
    while (true)
      if (parent->FirstChild()->ToComment() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
      else                                            break;
    
    //
    parent = XP_give_unique_expected_elem (parent, "VTKFile");
    stream->relink_downF();
    
    XP_check_expected_attribute (parent, "type", "UnstructuredGrid");
    XP_check_expected_attribute (parent, "version", "0.1");
    XP_check_expected_attribute (parent, "byte_order", "LittleEndian");
    
    //
    parent = XP_give_unique_expected_elem (parent, "UnstructuredGrid", false);
    stream->relink_downF();
    
    // read AppendedData section
    parent = parent->NextSibling();
    if (! parent)
      _errorr("No AppendedData section");
    
    if (_STRCMP(parent->Value(), "AppendedData") != 0)     _errorr ("name of xelement is not \"AppendedData\"");
    
    if (! (parent = parent->FirstChild ()))  _errorr("No child in AppendedData");
    if (   parent->Value()[0] != '_'      )  _errorr("There is no character \'_\' at the start of section AppendedData");
    parent = parent->NextSibling();
    
    if (_STRCMP(parent->Value(), "Input_Type") != 0)  _errorr("name of xelement is not Input_Type");

    int cll = this->file_type_s2e (parent->ToElement()->GetText());
    if (!add)
      this->set_coll (cll);
    else
      if (cll != this->coll)  errol;
    
    while (parent->NextSibling()) {
      errol;
    }
    
    // go to Piece
    XP_give_unique_expected_elem (stream->tixnod(), "Piece");
    stream->relink_downF();
    xnel = stream->tixel();
  }
  
  
  // variable declaration
  Stream auxstrm;
  // 
  bool status = true, prnt;
  int auxi;
  long i, j, nn, noel, cno, offset, auxl;
  char ed='-';
  const char *data=NULL, *offs=NULL, *typs, *str;
  char LINE[1023];
  int lLINE=1023;
  VTKrange range = VTKR_void;
  
  cno = 0;
  prnt = true;
  while (true) {
    //* key WORD
    if (lgc) {
      if (fgets (LINE, lLINE, stream->file()) == NULL)  break;
      str = LINE;
      SP_scan_word (str, WORD);
    }
    else {
      if (prnt) { if (! stream->relink_downF())  errol; else prnt = false; }
      else      { if (! stream->relink_next())   break; }
      WORD = (char *)stream->tixel()->Value();
    }
    
    //*  ******************
    //*  ***   POINTS   ***
    //*  ******************
    if      (_STRCMP(WORD, "Points") == 0) {
      if (range == VTKR_void)  range = VTKR_points;  else errol;
      
      // number of nodes
      if (lgc) {
    SP_scan_number (str, noel);
    SP_scan_expected_word2_exit  (str, "double", "float", "Invalid structure of the given VTK file, unsupported datatype of POINTS", CASE);
      }
      else {
    xnel->ToElement()->QueryLongAttribute("NumberOfPoints", &noel);
    data = XP_giveDAtext (stream->tixel(), 1, true, "ascii", "Float32", NULL, &(auxi=3));
      }
      
      if (!add) {
        // number of points
        this->set_and_alloc_nNodes (noel);

        // coords
        nNodes = 0;
        if (coll) { // 2d
          while (noel) { noel--;
            if (lgc) { status += Nodes[nNodes]->coords.scan_xyz (stream->file());  if (Nodes[nNodes]->coords[2] != 0.0)  _errorr("3. coordinates is not 0.0"); }
            else     { status += Nodes[nNodes]->coords.scan_xyz (data);            if (Nodes[nNodes]->coords[2] != 0.0)  _errorr("3. coordinates is not 0.0"); }
            nNodes ++;
          }
        }
        else {
          while (noel) { noel--;
            if (lgc) { status += Nodes[nNodes]->coords.scan_xyz (stream->file()); }
            else     { status += Nodes[nNodes]->coords.scan_xyz (data);            status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
            nNodes ++;
          }
        }
      }
      else {
        if (this->nNodes != noel)  errol;
        if (lgc) errol; // pokud legacy, tak tady musim preskocit radky asi noel radku
        else data = NULL;
      }
      
      if (lgc)  FP_skip_line (stream->file());
    }
    //*  *************************************
    //*  ***   UNSTRUCTURED GRID - CELLS   ***
    //*  *************************************
    else if (_STRCMP(WORD, "CELLS") == 0) {
      if (range == VTKR_points)  range = VTKR_cells;  else errol;
      
      // number of cells
      if (lgc) {
        SP_scan_number (str, noel);
        SP_scan_number (str, cno);
      }
      else {
        xnel->ToElement()->QueryLongAttribute("NumberOfCells" , &noel);
        data = XP_giveDAtext (stream->tixel(), 1, false, "ascii", "Int32", "connectivity", NULL);
        offs = XP_giveDAtext (stream->tixel(), 2, false, "ascii", "Int32", "offsets", NULL);
        typs = XP_giveDAtext (stream->tixel(), 3, true,  "ascii", "UInt8", "types", NULL);
      }
      
      if (!add) {
        this->set_and_alloc_nElems (noel);

        if (!lgc) offset = 0;
        // elem nodes
        for (i=0; i<noel; i++) {
          // nn
          if (lgc) { fscanf (stream->file(),"%ld",&nn); cno = cno - nn - 1; }
          else     { nn = offset;  status += SP_scan_number (offs, offset);  nn = offset - nn; }

          // Only strictly 2d or 3d meshes are supported.
          // So, the element type can be determined from number of nodes.
          int et;
          if (coll) { // 2D
            if (nn == 3)  et = 5;
            else errol;
          }
          else { // 3D
            if (nn == 8)  et = 12;
            else errol;
          }

          Elems[i]->set_type(et);

          //
          for (j=0; j<nn; j++)
            if (lgc) status += ( fscanf (stream->file(), "%ld", Elems[i]->nodes+j) == 1 );
            else     status += SP_scan_number (data, Elems[i]->nodes[j]);

        }

        if (lgc) FP_skip_line (stream->file(), 1);

        // types of elems
        if (lgc) {
          fgets (LINE, lLINE, stream->file());  str=LINE;
          SP_scan_expected_word_exit   (str, "CELL_TYPES", "Invalid structure of the given VTK file", CASE);
          SP_scan_expected_number_exit (str,         noel, "Invalid structure of the given VTK file, the number following CELL_TYPES");
        }

        for (i=0; i<noel; i++) {
          if (lgc)  status += ( fscanf (stream->file(), "%ld", &auxl) == 1);
          else      status += SP_scan_number (typs, auxl);
          if (auxl != Elems[i]->give_type())  errol;
        }
      }
      else {
        if (this->nElems != noel)  errol;
        if (lgc) errol; // pokud legacy, tak tady musim preskocit radky asi noel radku
        else { data = NULL; offs = NULL; }
      }
      
      if (lgc) FP_skip_line (stream->file(), 1);
    }
    //*  *******************************
    //*  ***   POLYDATA - ELEMENTS   ***
    //*  *******************************
    else if (_STRCMP(WORD, "LINES") == 0  ||  _STRCMP(WORD, "POLYGONS" ) == 0  ||  _STRCMP(WORD, "POLYS") == 0  ||  _STRCMP(WORD, "Vertices") == 0  ||  _STRCMP(WORD, "Verts") == 0  ||  _STRCMP(WORD, "Strips") == 0) {
      _errorr2("Invalid structure of the given VTK file, the key word %s is not supported", WORD);
    }
    //*  **********************
    //*  ***   POINT_DATA   ***
    //*  **********************
    else if (_STRCMP(WORD, "POINT_DATA") == 0  ||  _STRCMP(WORD, "PointData") == 0) {
      if (range == VTKR_cells  ||  range == VTKR_fields)  range = VTKR_pd;  else errol;
      ed = 'p';
      // head
      if (lgc) { SP_scan_expected_number_exit (str, nNodes, "The number following POINT_DATA is not equal to number of inclusions"); }
      else     { prnt = true; }
    }
    //*  *********************
    //*  ***   CELL_DATA   ***
    //*  *********************
    else if (_STRCMP(WORD, "CELL_DATA") == 0  ||  _STRCMP(WORD, "CellData") == 0) {
      if (range == VTKR_cells || range == VTKR_data)  range = VTKR_cd;  else errol;
      ed = 'c';
      // head
      if (lgc) { SP_scan_expected_number_exit (str, nElems, "The number following CELL_DATA is not equal to number of inclusions"); }
      else     { prnt = true; }
    }
    //*  **************************
    //*  ***   DATA - SCALARS   ***
    //*  **************************
    else if (_STRCMP(WORD, "SCALARS") == 0  ||  _STRCMP(WORD, "VECTORS") == 0  ||  _STRCMP(WORD, "TENSORS") == 0  ||  _STRCMP(WORD, "DataArray") == 0  ||  _STRCMP(WORD, "FIELD") == 0) {
      if (range == VTKR_pd || range == VTKR_cd || range == VTKR_data)  range = VTKR_data;  else errol;
      
      // blok dat se nacita jen pro FEM sit, proto zde kontroluju zda to neni mumech reseni
      //if (mtr != MTR_fem  &&  mtr != MTR_femadd)  continue;
      
      char type = WORD[0]; // save type of data - Scalar, Vector, Tensor, DataArray
      // head
      if (lgc)   SP_scan_word (str, WORD);
      else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      
      if (type != 'F') {
    //*  *** DATA for CELLS ***
    if (ed == 'c') {
      // bacic
      if      (_STRCMP(WORD, "INCLUSION_ID")       == 0) { auxl = 1;  scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'i', WORD);  for (j=0; j<nElems; j++)  { if (ST_scan_number (&auxstrm, Elems[j]->region) == false) _errorr2 ("There is few numbers in section %s", WORD);  Elems[j]->region--; } }
      else if (_STRCMP(WORD, "MATERIALS")          == 0) { auxl = 1;  scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'i', WORD);  for (j=0; j<nElems; j++)    if (ST_scan_number (&auxstrm, auxl            ) == false) _errorr2 ("There is few numbers in section %s", WORD); }
      else if (_STRCMP(WORD, "strain_planestrain") == 0) {
        stream->tixel()->ToElement()->QueryLongAttribute("NumberOfComponents", &auxl);
            this->set_elem_rslt_flags (false, true, false, pid, lc, 1);
        scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);
        for (j=0; j<nElems; j++) {
              Elems[j]->allocate_fields(pid);
              if (ST_scan_array(&auxstrm, auxl, Elems[j]->strain[pid][lc]) == false) _errorr2 ("There is few numbers in section %s", WORD);
              // termitovo: zde predpokladam voightovu notaci, proto delim dvojkou
              // zatim to je jen pro 2d, pro 3d bych prvne musel zkontrolovat, zda je dodrzena FEEP notace poradi clenu ve 3d, ve 2d je to jedno
              if (auxl != 3) errol; // now only 2d
              Elems[j]->strain[pid][lc][2] /= 2.0;
              if (auxl != P[pid]->nstrain_one()) {
                if(auxl == 3 && P[pid]->nstrain_one()==6) {
                  Elems[j]->strain[pid][lc][3] = Elems[j]->strain[pid][lc][2];
                  Elems[j]->strain[pid][lc][2] = Elems[j]->strain[pid][lc][4] = Elems[j]->strain[pid][lc][5] = 0.0;
        }
        else errol;
          }
            }
      }
      else if (_STRCMP(WORD, "stress_planestrain") == 0) {
        stream->tixel()->ToElement()->QueryLongAttribute("NumberOfComponents", &auxl);
            this->set_elem_rslt_flags (false, false, true, pid, lc, 1);
        scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);
        for (j=0; j<nElems; j++) {
              Elems[j]->allocate_fields(pid);
              if (ST_scan_array(&auxstrm, auxl, Elems[j]->stress[pid][lc]) == false) _errorr2 ("There is few numbers in section %s", WORD);
              if (auxl != P[pid]->nstrain_one()) {
                if(auxl == 3 && P[pid]->nstrain_one()==6) {
                  Elems[j]->stress[pid][lc][3] = Elems[j]->stress[pid][lc][2];
                  Elems[j]->stress[pid][lc][2] = Elems[j]->stress[pid][lc][4] = Elems[j]->stress[pid][lc][5] = 0.0;
        }
        else errol;
          }
            }
      }
      else
        _errorr2 ("Invalid name of SCALARS/DataArray for CELL_DATA/PointData \"%s\"", WORD);
      
    }
    //*  *** DATA for POINTS ***
    else if (ed == 'p') {
      if (_STRCMP(WORD, "uknw_displacement") == 0) {
        stream->tixel()->ToElement()->QueryLongAttribute("NumberOfComponents", &auxl);
        bool thzero = false;
        double auxd;
        
            if (auxl != P[pid]->ndisplc_one()) {
              if (auxl != 3  &&  P[pid]->ndisplc_one() != 2)
        errol;
          else
        thzero = true;
        }
        
            this->set_node_rslt_flags (true, false, false, pid, lc, 1);
        scan_DATA_field_head (stream, &auxstrm, str, type, auxl, 'f', WORD);
        if (thzero) auxl = 2;
            for (j=0; j<nNodes; j++) {
              Nodes[j]->allocate_fields(pid);
              if (ST_scan_array(&auxstrm, auxl, Nodes[j]->displc[pid][lc]) == false) _errorr2 ("There is few numbers in section %s", WORD);
              if (thzero) {
                ST_scan_number(&auxstrm,auxd);
                if (auxd != 0.0)
          errol;
              }
        }
      }
      else if (_STRCMP(WORD, "uknw_rotation") == 0) {
        if (lgc)  _errorr("have to skip values");
        
      }
      else
        _errorr2 ("Invalid name of SCALARS/DataArray for POINT_DATA/PointData \"%s\"", WORD);
    }
    else  errol;
      }
      else errol;
      
      if (lgc) FP_skip_line (stream->file(), 1);
      else { if (stream->tixel()->NextSibling() == NULL)  stream->relink_up(); }
    }
    else if (_STRCMP(WORD, "") == 0) {;}
    else _errorr2 ("Invalid structure of the given VTK file, unexpected key word \"%s\"", WORD);
    
    if (lgc) { if (cno) _errorr2 ("cno(%ld) is not zero", cno); }
    else {
      if (data && data[0] != '\0')  _errorr2("CellData of name \"%s\" has too many numbers", WORD);
      if (offs && offs[0] != '\0')  _errorr("error");
    }
  }
  
  if (lgc) delete [] WORD;
  if (status == false)  _errorr("error");
  
  //* close
  stream->close();
  delete stream;
  
}

void Mesh :: generate_regularSphereMesh_2d (int n, const double *a)
{
  if (!this->is_unset())  _errorr("Mesh in the function give_unitSphereMesh_3d is not empty/unset");
  if (n < 3 || n > 20)
    _errorr1 ("Number of segments must be between 3 and 20!");
  
  //
  this->set_coll(3);
  this->mtA = MT_GENERAL;
  
  // ALOCATING NODES
  this->set_and_alloc_nNodes (n * 4 + 1); // n nodes for each quater of inclusion * 4 quaters + 1 center point
  
  // SETTING RADIAL COORDS FOR NODES
  // >> Nodes[i = 0] = center of inclusion, already set to 0.0 in constructor
  // >> r = 1.0, no need to multiply by 'r'
  double phi = ((90.0 / n) * PI) / 180.0; // angle between two adjacent nodes (in radians)
  for (long i = 1, j = i - 1; i < this->nNodes; i++, j++) {
    this->Nodes[i]->coords.x = cos (phi * j);
    this->Nodes[i]->coords.y = sin (phi * j);
    // z not used and already set to 0.0 in constructor
  }
  
  // ALOCATING ELEMENTS
  this->set_and_alloc_nElems (n * 4); // n elements for each quater of inclusion * 4 quaters
  
  // SETTING CORRECT NODES FOR EACH ELEMENT
  int startNode = 1, node = startNode;
  for (long i = 0; i < this->nElems; i++) {
    this->Elems[i]->set_type (5); // = VTK_TRIANGLE
    this->Elems[i]->nodes[0] = 0; // all elements are connected with the Node[0] in center of inclusion
    this->Elems[i]->nodes[1] = node++; // nodes[1] = node, then node is incremented by 1
    this->Elems[i]->nodes[2] = (i == (this->nElems - 1)) ? startNode : node; // if last element then need co connect with startNode on first element
  }
  
  return;
}

void Mesh :: generate_regularSphereMesh_3d (int n, const double *a)
{
  if (!this->is_unset())  _errorr("Mesh in the function give_unitSphereMesh_3d is not empty/unset");
  if (n < 3 || n > 20)
    _errorr1 ("Number of segments must be between 3 and 20!");
  
  //
  this->set_coll(0);
  this->mtA = MT_GENERAL;
  
  // ALOCATING NODES
  int nNodesOnLayer = (n * 4); // n nodes for quater circle * 4 quaters
  int nLayersOfNodes = (n * 2 - 1); // n layers of nodes for hemisphere * 2 hemispheres - 1 duplicit center layer (top and bottom single nodes are not counted as layer)
  
  this->set_and_alloc_nNodes (nNodesOnLayer * nLayersOfNodes + 3);  // + 3 nodes (top, bottom and center)
  
  // SETTING SPHERICAL COORDS FOR NODES
  // r = 1.0, no need to multiply by 'r'
  
  // Nodes[0] = center of inclusion, already set to 0.0 in constructor
  this->Nodes[1]->coords.x = 1.0; // top node
  this->Nodes[this->nNodes - 1]->coords.x = -1.0; // bottom node
  
  double *phi = new double[2*n];
  // uniform distribution of angles
  if (a == NULL) {
    double ph = ((90.0 / n) * PI) / 180.0;
    for (long i = 0; i < 2*n; i++)
      phi[i] = ph;
  }
  else {
    phi[0] = PI/2.0 / n;
    double sum = PI/2.0 * a[1]/a[0];
    double k;
    
    do {
      phi[0] = phi[0] * 0.5 * PI / sum;
      k = phi[0] / give_curvature (0.5 * PI - 0, a[0], a[1]);
      sum = phi[0];
      for (long i = 1; i < n; i++) {
        phi[0+i] = k * give_curvature (0.5 * PI - sum, a[0], a[1]);
        sum += phi[0+i];
      }
    } while (fabs(sum - 0.5 * PI) > 0.01*PI);
    
    for (long i = 0; i < n; i++)   phi[n+i] = phi[i];
    for (long i = 0; i < n; i++)   phi[n-1-i] = phi[n+i];
    
    sum = 0.;
  }

  int nodeID = 2;
  double theta = ((90.0 / n) * PI) / 180.0; // angle between positive x-axis and node (in radians)

  double PHI = 0.0;
  for (long i = 0; i < nLayersOfNodes; i++) {
    PHI += phi[i];
    for (long j = 0; j < nNodesOnLayer; j++) {
      this->Nodes[nodeID]->coords.y = sin (PHI) * cos (theta * j);
      this->Nodes[nodeID]->coords.z = sin (PHI) * sin (theta * j);
      this->Nodes[nodeID]->coords.x = cos (PHI);
      nodeID++;
    }
  }

  delete [] phi;

  // ALOCATING ELEMENTS
  int nElemsInTopOrBottomLayers = n * 4; // n elements for each quater * 4 quaters
  int nElemsInMiddleLayer = n * 2 * 4; // n segments for each quater * 2 elements for segment (two tetrahedrons on top of each other) * 4 quaters
  int nMiddleLayers = (n - 1) * 2; // (n layers for hemisphere - top/bottom layer) * 2 hemispheres

  this->set_and_alloc_nElems (nElemsInTopOrBottomLayers * 2 + nElemsInMiddleLayer * nMiddleLayers);

  // SETTING CORRECT NODES FOR EACH ELEMENT
  int elemID = 0;

  // top layer
  int startNode = 2;
  nodeID = startNode;
  for (long i = 0; i < nElemsInTopOrBottomLayers; i++) {
    this->Elems[elemID]->set_type (10); // = VTK_TETRA
    this->Elems[elemID]->nodes[0] = 1; // all elements are connected with top node
    this->Elems[elemID]->nodes[1] = nodeID++; // nodes[1] = nodeID, then nodeID is incremented by 1
    this->Elems[elemID]->nodes[2] = (i == nElemsInTopOrBottomLayers - 1) ? startNode : nodeID; // if last element of layer then need co connect with startNode on first element
    this->Elems[elemID]->nodes[3] = 0; // all elements are connected with center node
    elemID++;
  }

  // middle layers
  int nextLayerNodeID, nextLayerStartNode;
  for (long i = 0; i < nMiddleLayers; i++) {
    startNode = nNodesOnLayer * i + 2; // n nodes on layer * number of previous layers + 2 (center and top)
    nodeID = startNode;
    nextLayerStartNode = startNode + nNodesOnLayer; // start node + n nodes on one layer
    nextLayerNodeID = nextLayerStartNode;
    for (long j = 0; j < nElemsInMiddleLayer; j += 2) { // step is two because two elems are set at the same time
      if (j % 4 == 0) { // first and all odd segments (to create zigzag pattern)
        this->Elems[elemID]->set_type (10); // = VTK_TETRA
        this->Elems[elemID]->nodes[0] = nodeID++;
        this->Elems[elemID]->nodes[1] = nextLayerNodeID;
        this->Elems[elemID]->nodes[2] = nodeID;
        this->Elems[elemID]->nodes[3] = 0; // center node
        elemID++;

        this->Elems[elemID]->set_type (10); // = VTK_TETRA
        this->Elems[elemID]->nodes[0] = nodeID;
        this->Elems[elemID]->nodes[1] = nextLayerNodeID++;
        this->Elems[elemID]->nodes[2] = nextLayerNodeID;
        this->Elems[elemID]->nodes[3] = 0; // center node
        elemID++;
      }
      else { // second and all even segments (to create zigzag pattern)
        this->Elems[elemID]->set_type (10); // = VTK_TETRA
        this->Elems[elemID]->nodes[0] = nodeID;
        this->Elems[elemID]->nodes[1] = nextLayerNodeID++;
        this->Elems[elemID]->nodes[2] = (j == (nElemsInMiddleLayer - 2)) ? nextLayerStartNode : nextLayerNodeID; // if last segment of layer then need co connect with nextLayerstartNode
        this->Elems[elemID]->nodes[3] = 0; // center node
        elemID++;

        this->Elems[elemID]->set_type (10); // = VTK_TETRA
        this->Elems[elemID]->nodes[0] = nodeID++;
        this->Elems[elemID]->nodes[1] = (j == (nElemsInMiddleLayer - 2)) ? nextLayerStartNode : nextLayerNodeID; // if last segment of layer then need co connect with nextLayerstartNode
        this->Elems[elemID]->nodes[2] = (j == (nElemsInMiddleLayer - 2)) ? startNode : nodeID; // if last segment of layer then need co connect with startNode
        this->Elems[elemID]->nodes[3] = 0; // center node
        elemID++;
      }
    }
  }

  // bottom layer
  startNode = this->nNodes - 1 - nNodesOnLayer; // all nodes - 1 (human vs computer numbering) - n nodes on one layer
  nodeID = startNode;
  for (long i = 0; i < nNodesOnLayer; i++) {
    this->Elems[elemID]->set_type (10); // = VTK_TETRA
    this->Elems[elemID]->nodes[0] = nodeID++;
    this->Elems[elemID]->nodes[1] = this->nNodes - 1; // last node
    this->Elems[elemID]->nodes[2] = (i == nNodesOnLayer - 1) ? startNode : nodeID;
    this->Elems[elemID]->nodes[3] = 0; // center node
    elemID++;
  }

  return;
}



//
void Mesh :: print_geometry_file_vtk (const char *rsltFileName, int pid, int lc, int nlc)
{
  if (nlc == 0) nlc = this->nLC;
  
  // Recursive call this function for separate load cases.
  if (nlc != 1) {
    // Loop over all load cases.
    for (int i=lc; i<lc+nlc; i++)
      this->print_geometry_file_vtk (rsltFileName, pid, i, 1);
    
    return;
  }
  
  
  long i, j, auxl;
  XMLElement *da;
  Stream *stream;
  bool lgc = true;
  bool twodim = P[pid]->give_twodim();
  char auxs[20001];  auxs[20000] = '\0'; 
  
  // allocate stream
  if (lgc)  stream = new Stream(STRM_file);
  else      stream = new Stream(STRM_tixel);
  
  //* START
  char *filename = new char[255];
  sprintf (filename, "%s.%02d.vtk", rsltFileName, lc); // do not delete filename, it is deleted inside of the function
  
  //
  print_VTK_START (stream, filename);
  
  //* HEAD
  char type[7];
  if (twodim) sprintf (type,"2Drec");
  else        sprintf (type,"3Drec");
  
  if (lgc)
    print_VTK_head (stream->file(), type);
  
  //* NODES
  da = print_VTK_nodes_head (stream, nNodes);
  
  for (i=0; i<nNodes; i++) {
    sprintf (auxs, "%15.8e %15.8e %15.8e", Nodes[i]->coords[0], Nodes[i]->coords[1], Nodes[i]->coords[2]);
    
    if (lgc) fprintf (stream->file(), "%s\n", auxs);
    else     da->InsertEndChild( stream->tix_doc()->NewText(auxs) );
  }
  
  // ELEMENTS
  if (lgc) {
    auxl = nElems;
    for (i=0; i<nElems; i++)  auxl += Elems[i]->nnodes;
    
    fprintf (stream->file(), "CELLS %ld %ld\n", nElems, auxl);
    for (i=0; i<nElems; i++) {
      fprintf (stream->file(), "%ld", Elems[i]->nnodes);
      for (j=0; j<Elems[i]->nnodes; j++)
    fprintf (stream->file(), " %ld", Elems[i]->nodes[j]);
      
      fprintf (stream->file(), "\n");
    }
    
    fprintf (stream->file(), "CELL_TYPES %ld\n", nElems);
    for (i=0; i<nElems; i++)
      fprintf (stream->file(), "%ld\n", Elems[i]->give_type());
    
  }
  else {
    _errorr("not implemented");
    // stream->tixel()->SetAttribute("NumberOfCells",  (int)Elems());
    // 
    // XMLElement *cells = stream->tix_doc()->NewElement("Cells");   stream->tixel()->InsertEndChild(cells);
    // 
    // XMLElement *dac = stream->tix_doc()->NewElement("DataArray");  cells->InsertEndChild(dac);
    // dac->SetAttribute("format", "ascii");
    // dac->SetAttribute("type", "Int32");
    // dac->SetAttribute("Name", "connectivity");
    // 
    // XMLElement *dao = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dao);
    // dao->SetAttribute("Name", "offsets");
    // 
    // XMLElement *dat = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dat);
    // dat->SetAttribute("type", "UInt8");
    // dat->SetAttribute("Name", "types");
    // 
    // char auxs[255];
    // long off=0;
    // for (long i=0; i<Elems(); i++) {
    //   Elems[i]->print_row_VTX (auxs);                                         dac->InsertEndChild(stream->tix_doc()->NewText(auxs));
    //   sprintf (auxs," %ld", off += Elems[i]->give_nno());                     dao->InsertEndChild(stream->tix_doc()->NewText(auxs));
    //   sprintf (auxs," %d", CellGeometry_e2i_VTK(Elems[i]->give_cellGeom()));  dat->InsertEndChild(stream->tix_doc()->NewText(auxs));
    // }
  }
  
  // POINT DATA
  print_VTK_init_point_data (stream, this->nNodes);
  
  if (no_dspl_rf && no_dspl_rf[pid] && no_dspl_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Displacemet", 'v', 'f', 0);    for (i=0; i<this->nNodes; i++) { print_values_Vector (auxs, this->Nodes[i]->displc[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  if (no_strn_rf && no_strn_rf[pid] && no_strn_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Strain",      't', 'f', 0);    for (i=0; i<this->nNodes; i++) { print_values_Tensor (auxs, this->Nodes[i]->strain[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  if (no_strs_rf && no_strs_rf[pid] && no_strs_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Stress",      't', 'f', 0);    for (i=0; i<this->nNodes; i++) { print_values_Tensor (auxs, this->Nodes[i]->stress[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  
  
  // CELL DATA
  print_VTK_init_cell_data (stream, this->nElems);
  
  if (el_dspl_rf && el_dspl_rf[pid] && el_dspl_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Displacemet", 'v', 'f', 0);    for (i=0; i<this->nElems; i++) { print_values_Vector (auxs, this->Elems[i]->displc[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  if (el_strn_rf && el_strn_rf[pid] && el_strn_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Strain",      't', 'f', 0);    for (i=0; i<this->nElems; i++) { print_values_Tensor (auxs, this->Elems[i]->strain[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  if (el_strs_rf && el_strs_rf[pid] && el_strs_rf[pid][lc]) {  da = print_VTK_data_head (stream, "Stress",      't', 'f', 0);    for (i=0; i<this->nElems; i++) { print_values_Tensor (auxs, this->Elems[i]->stress[pid][lc], twodim, 5);  print_auxs (lgc, stream, da, auxs); } }
  
  // termitovo tento if predelat na flag, a vubec tento radek je hrozne zbastlenej, udelat print_value_Scalar ....
  if (Elems[0]->error  &&  Elems[0]->error[pid]) { auxl = 1;   da = print_VTK_data_head (stream, "Error_elem",  's', 'f', auxl); for (i=0; i<this->nElems; i++) { print_values_Scalar (auxs, this->Elems[i]->error[pid]+lc,   auxl,   5);  print_auxs (lgc, stream, da, auxs); } }
  
  
  // END
  print_VTK_FINISH (stream);
  delete stream;
  
}//end of function: ()

void Mesh :: set_node_rslt_flags (bool displc_flag, bool strain_flag, bool stress_flag, int pid, int lc, int nlc)
{
  if (this->geom_is_readed() == false)  _errorr("problem geometry data not readed");
  if (nlc == 0)  _errorr("nlc == 0");
  
  if (displc_flag) { if (!no_dspl_rf)  no_dspl_rf = Allocate1Dbp(this->npa);  if (!no_dspl_rf[pid])  no_dspl_rf[pid] = Allocate1Dbz(nLC); }
  if (strain_flag) { if (!no_strn_rf)  no_strn_rf = Allocate1Dbp(this->npa);  if (!no_strn_rf[pid])  no_strn_rf[pid] = Allocate1Dbz(nLC); }
  if (stress_flag) { if (!no_strs_rf)  no_strs_rf = Allocate1Dbp(this->npa);  if (!no_strs_rf[pid])  no_strs_rf[pid] = Allocate1Dbz(nLC); }
  
  for (int i=lc; i<lc+nlc; i++) {
    if (displc_flag) { if (no_dspl_rf[pid][i] == false)  no_dspl_rf[pid][i] = true;  else _errorr("Mesh: displacement at nodes already computed for given load case");  }
    if (strain_flag) { if (no_strn_rf[pid][i] == false)  no_strn_rf[pid][i] = true;  else _errorr("Mesh: strain       at nodes already computed for given load case");  }
    if (stress_flag) { if (no_strs_rf[pid][i] == false)  no_strs_rf[pid][i] = true;  else _errorr("Mesh: stress       at nodes already computed for given load case");  }
  }
}
void Mesh :: set_elem_rslt_flags (bool displc_flag, bool strain_flag, bool stress_flag, int pid, int lc, int nlc)
{
  if (this->geom_is_readed() == false)  _errorr("problem geometry data not readed");
  if (nlc == 0)  _errorr("nlc == 0");

  if (displc_flag) { if (!el_dspl_rf)  el_dspl_rf = Allocate1Dbp(this->npa);  if (!el_dspl_rf[pid])  el_dspl_rf[pid] = Allocate1Dbz(nLC); }
  if (strain_flag) { if (!el_strn_rf)  el_strn_rf = Allocate1Dbp(this->npa);  if (!el_strn_rf[pid])  el_strn_rf[pid] = Allocate1Dbz(nLC); }
  if (stress_flag) { if (!el_strs_rf)  el_strs_rf = Allocate1Dbp(this->npa);  if (!el_strs_rf[pid])  el_strs_rf[pid] = Allocate1Dbz(nLC); }
  
  for (int i=lc; i<lc+nlc; i++) {
    if (displc_flag) { if (el_dspl_rf[pid][i] == false)  el_dspl_rf[pid][i] = true;  else _errorr("Mesh: displacement at elems already computed for given load case");  }
    if (strain_flag) { if (el_strn_rf[pid][i] == false)  el_strn_rf[pid][i] = true;  else _errorr("Mesh: strain       at elems already computed for given load case");  }
    if (stress_flag) { if (el_strs_rf[pid][i] == false)  el_strs_rf[pid][i] = true;  else _errorr("Mesh: stress       at elems already computed for given load case");  }
  }
}

void Mesh :: compute_node_fields (char flag, bool displc_flag, bool strain_flag, bool stress_flag, int pid, int lc, int nlc, PFCmode pfcMode)
{
  const Problem *prob = dynamic_cast <const Problem*>(P[pid]);
  if (prob == NULL) errol;
  
  nlc = Problem::check_lc_nlc(lc, nlc, nLC);
  this->set_node_rslt_flags (displc_flag, strain_flag, stress_flag, pid, lc, nlc);
  
  // this->allocare_fields_at_nodes
  for (int i=0; i<nNodes; i++)
    Nodes[i]->allocate_fields(pid);
  
  double coords[3];
  double **displc, **strain, **stress;
  
  for (int i=0; i<nNodes; i++) {
    Nodes[i]->coords.copy_to(coords);
    
    displc =  displc_flag ? Nodes[i]->displc[pid]+lc : NULL;
    strain =  strain_flag ? Nodes[i]->strain[pid]+lc : NULL;
    stress =  stress_flag ? Nodes[i]->stress[pid]+lc : NULL;
    
    prob->giveFieldsOfPoint(displc, strain, stress, coords, flag, lc, nlc, pfcMode, -3 , STRN_THEORETICAL_FEEP);
  }
}

void Mesh :: compute_element_fields (char flag, bool displc_flag, bool strain_flag, bool stress_flag, int pid, int lc, int nlc, PFCmode pfcMode)
{
  const Problem *prob = dynamic_cast <const Problem*>(P[pid]);
  if (prob == NULL) errol;
  
  nlc = Problem::check_lc_nlc(lc, nlc, nLC);
  this->set_elem_rslt_flags (displc_flag, strain_flag, stress_flag, pid, lc, nlc);
  
  //
  for (int i=0; i<nElems; i++)
    Elems[i]->allocate_fields(pid);
  
  for (int i=0; i<nElems; i++)
    Elems[i]->compute_centroids();
  
  double **displc, **strain, **stress;
  
  for (int i=0; i<nElems; i++) {
    displc =  displc_flag ? Elems[i]->displc[pid]+lc : NULL;
    strain =  strain_flag ? Elems[i]->strain[pid]+lc : NULL;
    stress =  stress_flag ? Elems[i]->stress[pid]+lc : NULL;
    
    Elems[i]->region = prob->giveFieldsOfPoint(displc, strain, stress, Elems[i]->centroids, flag, lc, nlc, pfcMode, Elems[i]->region, STRN_THEORETICAL_FEEP);
  }
}

double Mesh :: give_total_volume (void) const
{
  if (volume < 0.0) {
    double vol = 0.0;
    
    switch (this-> mtA) {
    case MT_REGULAR:
      if (this->coll)  vol = (point2[0]-point1[0])*(point2[1]-point1[1]);
      else              vol = (point2[0]-point1[0])*(point2[1]-point1[1])*(point2[2]-point1[2]);
      break;
    case MT_GENERAL:
      for (long i=0; i<nElems; i++)
    vol += this->Elems[i]->give_volume();
      break;
    default: _errorr( "Unknown mesh type");
    }
    
    ((Mesh*)this)->volume = vol;
  }
  
  return volume;
}


double Mesh :: give_regular_element_volume (void) const
{
  switch (this-> mtA) {
  case MT_REGULAR:  return  this->give_total_volume() / nElems;    break;
  case MT_GENERAL:  return  0.0;                                   break;
  default: _errorr( "Unknown mesh type");
  }
  
  return 0.0;
}

// ///
// void Mesh :: average_fields (bool displc_flag, bool strain_flag, bool stress_flag, int lc, int nlc)
// {
//   if (displc == NULL) { displc = new double*[nLC];  memset (displc, 0, nLC*sizeof(double*)); }
//   if (strain == NULL) { strain = new double*[nLC];  memset (strain, 0, nLC*sizeof(double*)); }
//   if (stress == NULL) { stress = new double*[nLC];  memset (stress, 0, nLC*sizeof(double*)); }
//   
//   double vol;
//   int ndisplc = P[0]->give_VECT_RANGE();
//   int nstrain = P[0]->give_VM_TENS_RANGE();
//   int nstress = P[0]->give_VM_TENS_RANGE();
//   
//   for (int i=lc; i<lc+nlc; i++) {
//     if (displc_flag  &&  elem_displc_rslt_flag[i] == false)  _errorr2("displc for lc %ld is not computed", i+1);
//     if (strain_flag  &&  elem_strain_rslt_flag[i] == false)  _errorr2("strain for lc %ld is not computed", i+1);
//     if (stress_flag  &&  elem_stress_rslt_flag[i] == false)  _errorr2("stress for lc %ld is not computed", i+1);
//     
//     if (displc_flag  &&  this->displc[i] == NULL)  displc[i] = new double[ndisplc];
//     if (strain_flag  &&  this->strain[i] == NULL)  strain[i] = new double[nstrain];
//     if (stress_flag  &&  this->stress[i] == NULL)  stress[i] = new double[nstress];
//     
//     if (displc_flag)  memset (displc[i], 0, ndisplc * sizeof(double));
//     if (strain_flag)  memset (strain[i], 0, nstrain * sizeof(double));
//     if (stress_flag)  memset (stress[i], 0, nstress * sizeof(double));
//     
//     for (long j=0; j<this->nElems; j++) {
//       vol = this->Elems[j]->give_volume();
//       
//       if (displc_flag)  AddVectorMultipledBy (vol, this->Elems[j]->displc[i], displc[i], ndisplc);
//       if (strain_flag)  AddVectorMultipledBy (vol, this->Elems[j]->strain[i], strain[i], nstrain);
//       if (stress_flag)  AddVectorMultipledBy (vol, this->Elems[j]->stress[i], stress[i], nstress);
//     }
//   }
// }


double Mesh :: give_homog_energy (int pid, int lc)
{
  const Problem *p = dynamic_cast <const Problem*>(P[pid]);
  if (p == NULL) errol;
  
  const double *hstrain = p->give_matrix()->give_globHomog_Strain(lc);
  const double *hstress = p->give_matrix()->give_globHomog_Stress(lc);
  
  
  double E;
  bool twodim = p->give_twodim();
  
  if (this->md == p->give_dimension()) {
    if (twodim)  E = MatrixOperations::giveTTproduct_1is2x2to3and2x2to3_SS (hstress, hstrain);
    else errol;
  }
  //        mesh is 2d           => planestrain is supposed
  else if ((this->coll == 3  &&  this->P[pid]->give_dimension() == 3)) {
    errol;
    //double strain2d[3], stress2d[3];
    //copy3Dto2Dtensors_FEEPreduced (hstrain, strain2d);
    //copy3Dto2Dtensors_FEEPreduced (hstress, stress2d);
    //
    //E = give_VectorVectorProduct (strain2d, stress2d, 3);
  }
  else errol;
  
  return E * this->give_total_volume();
}

double Mesh :: give_total_energy (int pid, int lc)
{
  if (el_strn_rf[pid][lc] == false  ||  el_strs_rf[pid][lc] == false)
    _errorr2("strain and stress for lc %ld is not computed", lc+1);
  
  double E = 0.0;
  
  for (long i=0; i<nElems; i++)
    E += Elems[i]->give_energy(pid, lc);
  
  return E;
}


double Mesh :: give_error (int pid1, int pid2, int lc)
{
  //double e2eA,e2eM, e2o;
  double strain11[3], stress11[3];
  double strain22[3], stress22[3];
  
  double vol, Eei,Eei2, Ee, Eai, Ea, Err_i;
  Ee = Ea = 0.0;
  
  for (long i=0; i<this->nElems; i++) {
    //    e2eA = thisAS->Elems[i]->energy[lc];
    //    e2eM = thisMM->Elems[i]->energy[lc];
    
    
    vol = this->Elems[i]->give_volume();
    
    this->Elems[i]->give_strain(strain11, pid1,lc);
    this->Elems[i]->give_strain(strain22, pid2,lc);
    
    this->Elems[i]->give_stress(stress11, pid1,lc);
    this->Elems[i]->give_stress(stress22, pid2,lc);
    
    // DEBUG
    // // stressM a stressA by mel byt stejny jako stresMM a stressAS
    // double stressM[3], stressA[3];
    // MatrixOperations::giveTVproduct_3is3x3to5and3(stressM, C, strainMM);
    // MatrixOperations::giveTVproduct_3is3x3to5and3(stressA, C, strainAS);
    // DEBUG
    
    SubtractVector(strain22, strain11, 3);
    SubtractVector(stress22, stress11, 3);
    
    // Energy error
    Eei = vol * give_VectorVectorProduct(strain11, stress11, 3);
    Ee += Eei;
    
    // Energy whole
    Eai = vol * give_VectorVectorProduct(strain22, stress22, 3);
    Ea += Eai;
    
    //
    Err_i = 100 * sqrt(Eei / Eai);
    this->Elems[i]->save_error(Err_i, 0, lc);
    

    // DEBUG
    if (Eai < 0.0) errol;    
    if (Eei/Eai < -1.0e-2) _errorr2("Eei = %e < 0.0", Eei/Eai);
    
    double C[5];
    this->Elems[i]->giveReducedStiffMatrix(C, pid1);
    Eei2 = vol * MatrixOperations :: giveVTVproduct_1is3and3x3to5and3(C, strain11);
    if (Eei2 < 0.0) errol;
    
    //if (abs((Eei2-Eei)/Eei) > 0.0001) {
    //  _warningg2("Eei  = %e < 0.0", Eei);
    //  _warningg2("Eei2 = %e < 0.0", Eei2);
    //  _warningg2("Eai = %e < 0.0",  Eai);
    //  //errol;
    //}
      
    // DEBUG
  }
  
  return 100 * sqrt(Ee / Ea);
}



} // end of namespace mumech

/*end of file*/