cell.cpp

Go to the documentation of this file.

#include "cell.h"

#include "point.h"
#include "attribute.h"

#include "gelib.h"
#include "arrays.h"
#include "mathlib.h"
#include "geometrylib.h"

#include <math.h>


namespace midaspace {

ComponentGeometry* allocComponentGeometry (classID cid, Cell *owner, ComponentGeometry* src=NULL) {
  switch (cid) {
  case classComponentGeometry1D           :  return (src ? new ComponentGeometry1D            (dynamic_cast<ComponentGeometry1D           *>(src)) : new ComponentGeometry1D           (owner, ZERO));  break;
  case classComponentGeometry2Dtriangle   :  return (src ? new ComponentGeometry2Dtriangle    (dynamic_cast<ComponentGeometry2Dtriangle   *>(src)) : new ComponentGeometry2Dtriangle   (owner, ZERO));  break;
  case classComponentGeometry2Dquadrangle :  return (src ? new ComponentGeometry2Dquadrangle  (dynamic_cast<ComponentGeometry2Dquadrangle *>(src)) : new ComponentGeometry2Dquadrangle (owner, ZERO));  break;
  case classComponentGeometry3Dhexahedron :  return (src ? new ComponentGeometry3Dhexahedron  (dynamic_cast<ComponentGeometry3Dhexahedron *>(src)) : new ComponentGeometry3Dhexahedron (owner, ZERO));  break;
  case classComponentGeometry3Dtetrahedron:  return (src ? new ComponentGeometry3Dtetrahedron (dynamic_cast<ComponentGeometry3Dtetrahedron*>(src)) : new ComponentGeometry3Dtetrahedron(owner, ZERO));  break;
  case classComponentGeometry2Dpolygon    :  return (src ? new ComponentGeometry2Dpolygon     (dynamic_cast<ComponentGeometry2Dpolygon    *>(src)) : new ComponentGeometry2Dpolygon    (owner, ZERO));  break;
  case classComponentGeometry1Dpoly       :  return (src ? new ComponentGeometry1Dpoly        (dynamic_cast<ComponentGeometry1Dpoly       *>(src)) : new ComponentGeometry1Dpoly       (owner, ZERO));  break;
  default: errol;
  }
  return NULL;
}

//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS CELL                          *** *** *** ***
//*  ********************************************************************************************
Cell :: Cell (classID mecg, long gid, long oid, long prop, const Geometry *mg, long nn, long ne, long nf) : GeometryComponent (mg, gid, oid, prop) //, ComponentGeometry()
{
  // this has to be before cg
  ;       points.resizeup(nn);
  if (ne) edges. resizeup(ne);
  if (nf) faces. resizeup(nf);
  
  if (mecg)  cg = allocComponentGeometry (mecg, this);
  else       cg = NULL;
  
  //
  connectivity_assembled = false;
  //
  duplmaster = NULL;
}
Cell :: Cell (classID mecg, const Cell * src) : GeometryComponent (src)
{
  //
  points.be_copy_of(&src->points);
  edges.resizeup(src->edges());
  faces.resizeup(src->faces());
  //
  connectivity_assembled = false;
  //
  duplmaster = NULL;
  
  if (((Mesh*)Geom)->connectivity_is_assembled())  _errorr ("error");
  
  if (src->cg && src->cg->give_classid() == mecg) { cg = allocComponentGeometry (mecg, this, src->cg);  cg->set_owner(this); }
  else errol;
}
Cell :: ~Cell()
{
  delete cg;
}


//*  *********************************************
//*  ***   ***   ***   DUPLICITY   ***   ***   ***
//*  *********************************************
void Cell :: assure_duplicity_master (void)
{
  if ( duplmaster != (Cell *)-1 )  return;
  if ( dynamic_cast<FElement *>(this) == NULL )  _errorr ("Error");
  check_connectivity();

  int i;
  FElement *master;

  for (i=0; i<edges[0]->give_numsuperelem(); i++) {
    master = (FElement*)edges[0]->give_superelem(i);
    if (master == this) continue;

    if (master->has_same_geom_with(this)) {
      duplmaster = master;
      if ( this->give_ID() < duplmaster->give_ID() )  _errorr ("Ascending duplicity chain");
      break;
    }
  }
  if (i == edges[0]->give_numsuperelem())
    duplmaster = NULL;

}

Cell* Cell :: last_duplicity_master (void) const
{
  if (this == duplmaster)
    _errorr ("");
  
  ((Cell*)this)->assure_duplicity_master();
  
  if (!duplmaster)  return NULL;

  Cell *master;
  if ( (master = duplmaster->last_duplicity_master()) )  return master;

  return duplmaster;
}

void Cell :: integrate_duplicated_cell (const Cell *slave)
{
# ifdef DEBUG
  //
  if (!points.has_same_members_as(slave->points))  _errorr ("Duplicity slave and master differ in list of nodes");
  long i;
  if (edges() != slave->edges()) errol;
  for (i=0; i<edges(); i++)
    if ( (       edges[i]->last_duplicity_master() ?        edges[i]->last_duplicity_master() :        edges[i]) !=
     (slave->edges[i]->last_duplicity_master() ? slave->edges[i]->last_duplicity_master() : slave->edges[i])    )
      _errorr ("Duplicity slave and master differ in list of edges");
  if (faces() != slave->faces()) errol;
  for (i=0; i<faces(); i++)
    if ( (       faces[i]->last_duplicity_master() ?        faces[i]->last_duplicity_master() :        faces[i]) !=
     (slave->faces[i]->last_duplicity_master() ? slave->faces[i]->last_duplicity_master() : slave->faces[i])    )
      _errorr ("Duplicity slave and master differ in list of faces");
# endif
}

void Cell :: switch_node_pointer (Point *slave, Point *master, bool duplcheck)
{
  if (points.replace_member_by (slave, master))
    if ( this->check_collapse() )
      duplmaster = this;
}
bool Cell :: switch_edge_pointer (Edge *slave, Edge *master)
{
  return edges.replace_member_by (slave, master);
}
bool Cell :: switch_face_pointer (Face *slave, Face *master)
{
  return faces.replace_member_by (slave, master);
}


void Cell :: setup_duplicity_master (Cell *newmaster)
{
  if ( newmaster == ((Cell*)-1) ) {
    if (duplmaster)  _errorr ("");
    duplmaster = newmaster;
    return;
  }

  if (this == duplmaster)
    errol;

  if ( this->give_classid() != newmaster->give_classid() )
    _errorr ("Mismatch");

  if ( newmaster == duplmaster )  return;

  Cell *oldmaster = this->last_duplicity_master();

  if ( newmaster == oldmaster )  return;

  if (oldmaster == NULL) {
    duplmaster = newmaster;
    if ( this->give_ID() < duplmaster->give_ID() )  _errorr ("Ascending duplicity chain");
  }
  else
    if ( newmaster->give_ID() > oldmaster->give_ID() )
      newmaster->setup_duplicity_master(oldmaster);
    else
      oldmaster->setup_duplicity_master(newmaster);

}

bool Cell :: invisible_duplicated (char flag)
{
  if  (flag != 'a' && flag != 'c' && flag != 'r')  _errorr("forbidden flag");

  this->assure_duplicity_master();

  // no duplicity
  if (duplmaster == NULL)  return false;

  // collapsed cell
  if (duplmaster == this) {
    if (flag == 'r')  return false;
    this->connectivity_removing();
  }
  // duplicated
  else {
    if (flag == 'c')  return false;
    this->last_duplicity_master()->integrate_duplicated_cell(this);
    this->switch_myself_at_connectivity(this->last_duplicity_master());
  }
  
  return true;
}

void Cell :: switch_myself_at_connectivity (Cell *master)
{
# ifdef DEBUG
  if (!master) errol;
# endif
  this->connectivity_removing();
}


//*  ********************************************
//*  ***   ***   ***   ASSORTED   ***   ***   ***
//*  ********************************************

bool Cell :: has_same_geom_with (Cell *slave) const
{
  if (this == slave)   return false;
  //return nodes.has_same_members_as(slave->nodes);
  return points.has_these_members(slave->give_points()[0]);
}

bool Cell :: cross_abscissa_node (const PoinT *a1, const PoinT *a2, long cunn, const Point **unnod, const Point *&nod, double &ksi, PoinT *cp) const
{
  long i, j;
  for (i=0; i<points(); i++) {
    for (j=0; j<cunn; j++)
      if (points[i] == unnod[j]) break;

    if (j==cunn && points[i]->give_coords()->give_ksiAtAbscissa(ZERO, this->give_circum(), a1, a2, ksi)) {
      nod = points[i];
      cp->bePointAtAbscissa (a1, a2, ksi);
      return true;
    }
  }
  return false;
}

int Cell :: cross_abscissa_face (const PoinT *a1, const PoinT *a2, long cunf, const Face **unfac, const Cell *&comp, PoinT *nc, PoinT *cp)
{
  // napsano pro brick
  long i,j,nr;
  double ksi[2],eta[2],t[2];
  PoinT i1, i2, *ip;

  for (i=0; i<faces(); i++) {
    for (j=0; j<cunf; j++)  if (faces[i] == unfac[j]) break;
    if (j!=cunf) continue;

    nr =  intersec_rectangle3d_line (ZERO, this->give_circum(),
                     faces[i]->points[0]->give_coords(),faces[i]->points[1]->give_coords(),
                     faces[i]->points[2]->give_coords(),faces[i]->points[3]->give_coords(),
                     a1,a2,ksi,eta,t,&i1,&i2);
    if (nr==0)  continue;
    if (nr==-1)
      _errorr("Something is wrong in function cross_abscissa_face");

    if ( fabs(ksi[0])<(1.0+ZERO) && fabs(eta[0])<(1.0+ZERO) && fabs(t[0])<(1.0+ZERO) ){
      nc->x = ksi[0];
      nc->y = eta[0];
      ip    = &i1;
      if (nr==1) break;
      if (nr==2) nr=3;
    }

    if (nr==1) continue;

    if ( fabs(ksi[1])<(1.0+ZERO) && fabs(eta[1])<(1.0+ZERO) && fabs(t[1])<(1.0+ZERO) ){
      if (nr==3) continue;
      if (nr==2){
    nc->x = ksi[1];
    nc->y = eta[1];
    ip    = &i2;
    break;
      }
    }

    if (nr==3) break;
  }

  if (i==faces())  return 0;

  const Face *fac = faces[i];

  cp->copy(ip);

  const Edge *edg;
  const Point *t1, *t2;

  if      ( fabs(nc->x-1.0) < ZERO ) { t1=fac->give_point(3); t2=fac->give_point(0); nc->x=nc->y; edg=(Edge*)fac->give_edge(3); }  // ksi =  1.0
  else if ( fabs(nc->x+1.0) < ZERO ) { t1=fac->give_point(2); t2=fac->give_point(1); nc->x=nc->y; edg=(Edge*)fac->give_edge(1); }  // ksi = -1.0
  else if ( fabs(nc->y-1.0) < ZERO ) { t1=fac->give_point(1); t2=fac->give_point(0);              edg=(Edge*)fac->give_edge(0); }  // eta =  1.0
  else if ( fabs(nc->y+1.0) < ZERO ) { t1=fac->give_point(2); t2=fac->give_point(3);              edg=(Edge*)fac->give_edge(2); }  // eta = -1.0
  else {
    comp = fac;
    return 2;
  }

  comp = edg;
  ((Edge*)edg)->transform_nc (t1,t2,nc->x);

  return 1;
}



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS FACEDGE                         *** *** *** ***
//*  ********************************************************************************************
void Facedge :: setadd_masterel (const Element *comp)
{
  int dim = comp->give_dimension();
  
  if (dim > masterel_dim)  return;
  
  if (dim < masterel_dim) {
    masterel_dim = dim;
    masterels.resize_zero();
  }
  
  masterels.add(comp);
}

const Element* Facedge :: give_main_masterel_uniq (void) const
{
  check_connectivity();
  
  if (masterel_dim == this->give_dimension()) {
    if (masterels() == 1)
      return masterels[0];
    else
      _errorr3("Facedge %ld of dim %d has multiple superelems with same dimension", this->give_ID(), this->give_dimension());
  }
  
  return NULL;
}

void Facedge :: integrate_duplicated_cell (const Facedge *slave)
{
  Cell :: integrate_duplicated_cell (slave);
  
  superelems.add_unique (slave->give_superelems()[0]);
}




//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS EDGE                          *** *** *** ***
//*  ********************************************************************************************
void Edge :: connectivity_assembling (bool re)
{
  if (!re && this->connectivity_assembled == true) return;
  this->connectivity_assembled = true;
  
  int i;
  for (i=0; i<points(); i++) const_cast<Point*>(points[i])->setadd_superedge(this, re);
}
void Edge :: connectivity_removing (void)
{
  if (this->connectivity_assembled == false) return;
  this->connectivity_assembled = false;
  
  int i;
  for (i=0; i<points(); i++) const_cast<Point*>(points[i])->remove_superedge(this);
}

void Edge :: switch_myself_at_connectivity (Cell *master)
{
  Cell :: switch_myself_at_connectivity(master);
  
  for (int i=0; i<superfaces(); i++)  const_cast<Face*   >(superfaces[i])->switch_edge_pointer (this, (Edge*)master);
  for (int i=0; i<superelems(); i++)  const_cast<Element*>(superelems[i])->switch_edge_pointer (this, (Edge*)master);
}

void Edge :: set_model_prop (long val, const Model *model, bool flag)
{
  const Edge *mdl_master_edge;
  
  if (flag == false) {
    if (val > model->give_Nels())
      val = val - model->give_Nels();
    else
      flag = true;
  }
  
  if (flag) {
# ifdef DEBUG
    if (model->give_Elem(val-1)->give_ned() != 1) errol;
# endif
    mdl_master_edge = model->give_Elem(val-1)->give_edge(0);
    val = mdl_master_edge->give_ID()+1;
  }
  else {
    mdl_master_edge = model->give_Edge(val-1);
  }
  
  if ( this->checkset_mprop(val) )
    this->give_facedgeAttribs()->set_new (mdl_master_edge->give_facedgeAttribs());
  
}

void Edge :: integrate_duplicated_cell (const Edge *slave)
{
  Facedge :: integrate_duplicated_cell (slave);
  
  superfaces.add_unique (slave->give_superfaces()[0]);
}


/*
void edge::transform_nc (long t,double ksi)
// releas verse
{
  if (t==nodes[1]) ksi = -ksi;
}
*/
void Edge :: transform_nc (const Point *t1, const Point *t2, double &ksi)
// debug verse
{
  if      (t1==points[1] && t2==points[0])  ksi = -ksi;
  else if (t1!=points[0] || t2!=points[1])  exit (2);
}
bool Edge :: is_nod_on (double norm, const PoinT *point, double &ksi) const
{
  return point->give_ksiAtAbscissa (ZERO, norm, points[0]->give_coords(), points[1]->give_coords(), ksi);
}

const Face* Edge :: give_superface (long ne, const Edge **edgs) const
{
  for (int i=0; i<superfaces(); i++)
    if (superfaces[i]->give_ned() == ne)
      if ( superfaces[i]->give_edges()->has_these_members(ne, edgs) )
    return superfaces[i];

  return NULL;
}

void Edge :: switch_node_pointer (Point *slave, Point *master, bool duplcheck)
{
  Cell :: switch_node_pointer (slave, master, duplcheck);

  if (duplmaster == this  ||  !duplcheck)  return;

  for (int i=0; i<master->give_numsuperedge(); i++)
    if ( master->give_superedge(i)->has_same_geom_with(this) ) {
      if ( this->give_ID() > master->give_superedge(i)->give_ID() )
    this->setup_duplicity_master((Edge*)master->give_superedge(i));
      else
    ((Edge *)master->give_superedge(i))->setup_duplicity_master(this);
    }
}

void Edge :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  switch (fff) {
  case FFF_T3d: {
    const Element* master = this->give_main_masterel_uniq();
    
    // name + id
    fprintf (stream, "curve %3ld   vertex ", this->give_ID()+1);
    
    // vertices
    if (points() != 2)  errol;
    fprintf (stream, " %3ld %3ld", points[0]->give_ID()+1, points[1]->give_ID()+1);
    
    // property
    fprintf (stream, "   property %ld", master ? master->give_ID() + 1 : Geom->give_Nels() + this->give_ID() + 1);
    
    // size
    double size = this->give_elemSize();
    if (master && master->isTruss()) {
      if (superelems() !=1) errol;
      size = -1;
    }
    if (size > 0.0)  fprintf (stream, "   size %lf", size);
    else             fprintf (stream, "   count %d", (int)-size);
    
    if (master) {
      const GelemAttribs *gatt = dynamic_cast<const GelemAttribs*>(master->give_elemAttribs());
      if (!gatt)  errol;
      
      if (gatt->give_virtual())  fprintf (stream, "   output no");
      else                       fprintf (stream, "   output yes");
    }
    break;
  }
  default: errol;
  }
  
  //
  if (endline)  fprintf (stream, "\n");
}

double Edge :: give_elemSize (void) const
{
  double size;
  const Element* master = this->give_main_masterel_uniq();
  
  if (master)
    size = ((FacedgeAttribs*)this->attributes)->give_elemSize();
  else
    size = superfaces[0]->give_elemSize();
  
  if (master) {
    long count = ((FacedgeAttribs*)this->attributes)->give_elemCount();
    if (count) {
      double cntsize = master->give_characteristic_size() / count;
      if (size > cntsize)
    size = - count;
    }
  }
  
  if (master) {
    const GelemAttribs *gatt = dynamic_cast<const GelemAttribs*>(master->give_elemAttribs());
    if (!gatt)  errol;
    if (! gatt->give_virtual()) {
      if (gatt->give_mat()->give_type() == MAT_RIGID)
    size = -1;
    }
  }
  
  return size;
}


//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS FACE                          *** *** *** ***
//*  ********************************************************************************************
void Face :: connectivity_assembling (bool re)
{
  if (!re && this->connectivity_assembled == true) return;
  this->connectivity_assembled = true;
  
  int i;
  for (i=0; i<points(); i++) const_cast<Point*>(points[i])->setadd_superface(this, re);
  for (i=0; i<edges(); i++)  const_cast<Edge*> (edges[i]) ->setadd_superface(this, re);
}
void Face :: connectivity_removing (void)
{
  if (this->connectivity_assembled == false) return;
  this->connectivity_assembled = false;
  
  int i;
  for (i=0; i<points(); i++) const_cast<Point*>(points[i])->remove_superface(this);
  for (i=0; i<edges(); i++)  const_cast<Edge*> (edges[i]) ->remove_superface(this);
}

void Face :: switch_myself_at_connectivity (Cell *master)
{
  Cell :: switch_myself_at_connectivity(master);
  
  for (int i=0; i<superelems(); i++)  const_cast<Element*>(superelems[i])->switch_face_pointer (this, (Face*)master);
}

void Face :: set_model_prop (long val, const Model *model, bool flag)
{
  const Face *mdl_master_face;
  
  if (flag == false) {
    if (val > model->give_Nels())
      val = val - model->give_Nels();
    else
      flag = true;
  }
  
  if (flag) {
# ifdef DEBUG
    if (model->give_Elem(val-1)->give_nfa() != 1) errol;
# endif
    mdl_master_face = model->give_Elem(val-1)->give_face(0);
    val = mdl_master_face->give_ID()+1;
  }
  else {
    mdl_master_face = model->give_Face(val-1);
  }
  
  if ( this->checkset_mprop(val) )
    this->give_facedgeAttribs()->set_new (mdl_master_face->give_facedgeAttribs());
  
}


/*
void face::transform_nc (const long *t,double *nc)
// releas verse
{
  if      (t1==edges[0]) { nc->z=nc->y; nc->y= nc->x; if (t2==edges[1]) nc->x=-nc->z;
                                                        else                nc->x= nc->z; }
  else if (t1==edges[1]) {              nc->x=-nc->x; if (t2==edges[2]) nc->y=-nc->y; }
  else if (t1==edges[2]) { nc->z=nc->y; nc->y=-nc->x; if (t2==edges[3]) nc->x= nc->z;
                                                        else                nc->x=-nc->z; }
  else                                                  if (t2!=edges[0]) nc->y=-nc->y;
}
*/
void Face :: transform_nc (const Edge *t1, const Edge *t2, PoinT *nc) const
// debug verse
{
  if      (t1==edges[0]) { nc->z=nc->y; nc->y= nc->x; if      (t2==edges[1]) nc->x=-nc->z;
                                                      else if (t2==edges[3]) nc->x= nc->z;
                                                      else                   exit (3);     }
  else if (t1==edges[1]) {              nc->x=-nc->x; if      (t2==edges[2]) nc->y=-nc->y;
                                                      else if (t2!=edges[0]) exit (4);     }
  else if (t1==edges[2]) { nc->z=nc->y; nc->y=-nc->x; if      (t2==edges[3]) nc->x= nc->z;
                                                      else if (t2==edges[1]) nc->x=-nc->z;
                                                      else                   exit (5);     }
  else if (t1==edges[3]) if        (t2==edges[2]) nc->y=-nc->y;
                         else { if (t2!=edges[0]) exit (6); }
  else exit (7);
}


void Face :: switch_node_pointer (Point *slave, Point *master, bool duplcheck)
{
  Cell :: switch_node_pointer (slave, master, duplcheck);

  if (duplmaster == this  ||  !duplcheck)  return;

  for (int i=0; i<master->give_numsuperface(); i++)
    if ( master->give_superface(i)->has_same_geom_with(this) ) {
      if ( this->give_ID() > master->give_superface(i)->give_ID() )
    this->setup_duplicity_master((Face*)master->give_superface(i));
      else
    ((Face *)master->give_superface(i))->setup_duplicity_master(this);
    }
}

void Face :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  switch (fff) {
  case FFF_T3d: {
    const VectoR *norm = NULL;
    
    switch (cg->give_elemGeom()) {
    case CG_Triangle:
    case CG_Quadrangle:
    case CG_Polygon: {
      const Gelement* master = dynamic_cast<const Gelement*>(this->give_main_masterel_uniq());
      if (!master) errol;
      
      // name + id
      fprintf (stream, "patch %3ld", this->give_ID()+1);
      
      // normal
      norm = ((ComponentGeometry2D*)cg)->give_normal();
      fprintf (stream, "   normal %21.14e %21.14e %21.14e", norm->x, norm->y, norm->z);
      
      // property
      fprintf (stream, "   property %ld", master ? master->give_ID() + 1 : Geom->give_Nels() + this->give_ID() + 1);
      
      // boundary curves
      fprintf (stream, "   boundary curve ");
      int i=0, frwd=1;
      for (i=0; i<edges()-1; i++) {
    if      (edges[i]->give_point(0) == points[i  ]  &&  edges[i]->give_point(1) == points[i+1])  frwd =  1;
    else if (edges[i]->give_point(0) == points[i+1]  &&  edges[i]->give_point(1) == points[i  ])  frwd = -1;
    else    errol;
    fprintf (stream, " %3ld", frwd*(edges[i]->give_ID()+1));
      }
      if      (edges[i]->give_point(0) == points[i]  &&  edges[i]->give_point(1) == points[0])  frwd =  1;
      else if (edges[i]->give_point(0) == points[0]  &&  edges[i]->give_point(1) == points[i])  frwd = -1;
      else    errol;
      fprintf (stream, " %3ld", frwd*(edges[i]->give_ID()+1));
      
      // size
      fprintf (stream, "   size %lf", this->give_elemSize());
      
      // count
      if (master) {
        const GelemAttribs *gatt = dynamic_cast<const GelemAttribs*>(master->give_elemAttribs());
        if (!gatt)  errol;
        
        if (gatt->give_virtual())  fprintf (stream, "   output no");
        else                       fprintf (stream, "   output yes");
      }
      
      // fixed
      if (master) {
    const GPA<const Vertex> *fixvers = master->give_fixedVertices();
    if (fixvers) {
      fprintf (stream, "  fixed vertex");
      for (i=0; i<fixvers[0](); i++)
        fprintf (stream, " %ld", fixvers[0][i]->give_ID()+1);
    }
    const GPA<const Gelement> *fixgelems = master->give_fixedGelements();
    if (fixgelems) {
      fprintf (stream, "  fixed curve");
      for (i=0; i<fixgelems[0](); i++)
        fprintf (stream, " %ld", fixgelems[0][i]->give_same_dimension_facedge()->give_ID()+1);
    }
      }
      
      //
      if (Pd->give_PDBO(PDBO_P_mesh_quads)) 
    fprintf (stream, "   quad map yes equidistant");
      
      break;
    }
    default: errol;
    }
    break;
  }
  default: errol;
  }
  
  //
  if (endline)  fprintf (stream, "\n");
}

double Face :: give_elemSize (void) const
{
  double size;
  const Element* master = this->give_main_masterel_uniq();
  
  //* size
  if (master)
    size = ((FacedgeAttribs*)this->attributes)->give_elemSize();
  else
    errol;
    // tady bude asi volani volume
  
  //* Rigid Material
  if (master) {
    const GelemAttribs *gatt = dynamic_cast<const GelemAttribs*>(master->give_elemAttribs());
    if (!gatt)  errol;
    if (! gatt->give_virtual()) {
      if (gatt->give_mat() && gatt->give_mat()->give_type() == MAT_RIGID)
    size = master->give_characteristic_size();
    }
  }
  
  return size;
}


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS ELEMENT                         *** *** *** ***
//*  ********************************************************************************************
VTKPDtopology Element :: give_VTKPDtopology(void) const
{
  switch (CellGeometry_e2i_VTK(this->give_cellGeom())) {
  case 3:   return VTKPD_LINE;      // VTK_LINE
  case 5:   return VTKPD_POLYGON;   // VTK_TRIANGLE
  case 9:   return VTKPD_POLYGON;   // VTK_QUAD
  default:  _errorr ("element is not polydata unsuported VTK cell type");
  }
  return VTKPD_Void;
}
// ///
// const Facedge* Element :: give_same_dimension_facedge (void) const
// {
//   switch (this->give_dimension()) {
//   case 1:  return edges[0];
//   case 2:  return faces[0];
//   case 3:   errol;
//   default:  errol;
//   }
//   return NULL;
// }
const Facedge* Element :: give_same_dimension_facedge (void) const
{
  switch (this->give_dimension()) {
  case 1:  if (edges() != 1) errol;  return edges[0];
  case 2:  if (faces() != 1) errol;  return faces[0];
  case 3:   errol;
  default:  errol;
  }
  return NULL;
}


void Element :: connectivity_assembling (bool re)
{
  if (!re && this->connectivity_assembled == true) return;
  this->connectivity_assembled = true;
  
  //*  *** NODES ***
  for (int i=0; i<points(); i++)
    const_cast<Point*>(points[i])->setadd_superelem(this, re);
  
  long i, nned, nnfa;
  bool all_exist_subs;
  const Point **subnodes = NULL;  // nodes at sub-element
  const Edge  **subedges = NULL;  // edges at sub-element
  //
  Edge *edg;
  Face *fac;
  
  
  //*  *** EDGES ***
  nned = this->give_edge_nodes (subnodes);
  all_exist_subs = true;
  
  for (i=0; i<edges(); i++) {
    // edge uz je prirazeno
    if (edges[i]) {
      if (!re) errol; // to lze jen kdyz dam re
      if (!edges[i]->give_points()->has_these_members(nned, subnodes+nned*i) )  errol;   // kontrola ze uz prirazena edge ma spravny nodes
      edg = (Edge*)edges[i];
    }
    // priradim nejakou existujici, nebo nove naallocovanou
    else {
      edg = (Edge*) subnodes[nned*i]->give_superedge (nned, subnodes+nned*i);
      if (edg==NULL) {
    edg = new Edge(nned, subnodes+nned*i);
    const_cast<Geometry*>(Geom)->add_another_Edge (edg);
    all_exist_subs = false;
      }
      this->set_edge(i, edg);
    }
    
    edg->connectivity_assembling(re);
    
    edg->setadd_superelem(this, re);
  }
  if (all_exist_subs)  this->setup_duplicity_master ((Element*)-1);
  
  delete [] subnodes;  subnodes = NULL;
  
  
  nnfa = this->give_face_nodes_edges (subnodes, subedges);
  all_exist_subs = true;
  
  for (i=0; i<faces(); i++) {
    // face uz je prirazeno
    if (faces[i]) {
      if (!re) errol; // to lze jen kdyz dam re
      if (!faces[i]->give_edges()->has_these_members(nnfa, subedges+nnfa*i) )  errol;   // kontrola ze uz prirazena face ma spravny edges
      fac = (Face*)faces[i];
    }
    // priradim nejakou existujici, nebo nove naallocovanou
    else {
      fac = (Face*) subedges[nnfa*i]->give_superface (nnfa, subedges+nnfa*i);
      if (fac==NULL) {
    fac = new Face(nnfa, nnfa, subnodes+nnfa*i, subedges+nnfa*i);
    const_cast<Geometry*>(Geom)->add_another_Face (fac);
    all_exist_subs = false;
      }
      this->set_face(i, fac);
    }
    
    fac->connectivity_assembling(re);
    
    fac->setadd_superelem(this, re);
  }
  if (nnfa && all_exist_subs)  if (duplmaster == NULL) errol;
  
  delete [] subnodes;
  delete [] subedges;
}
void Element :: connectivity_removing (void)
{
  if (this->connectivity_assembled == false) return;
  this->connectivity_assembled = false;
  
  int i;
  for (i=0; i<points(); i++) const_cast<Point*>(points[i])->remove_superelem (this);
  for (i=0; i<edges(); i++)  const_cast<Edge*> (edges[i]) ->remove_superelem (this);
  for (i=0; i<faces(); i++)  const_cast<Face*> (faces[i]) ->remove_superelem (this);
  
  // toto smaz pokut to zde nebude padat
  // for (i=0; i<points(); i++)                const_cast<Point*>(points[i])->remove_superelem (this);
  // for (i=0; i<edges(); i++)  if (edges[i])  const_cast<Edge*> (edges[i]) ->remove_superelem (this);  else if (master) errol;
  // for (i=0; i<faces(); i++)  if (faces[i])  const_cast<Face*> (faces[i]) ->remove_superelem (this);  else if (master) errol;
}

void Element :: set_load (int i, int indx)
{
  const BoundaryCond *bc = Pd->give_BC(i);  if (bc == NULL) _errorr2("There is no boundary condition with ID = %d", i+1);
  
  switch (bc->give_loctype()) {
  case BC_CEL:  this->connectivity_assembling(false);  ((Geometry*)Geom)->assign_EdgeLoad_by_elem_ID(bc, this->give_ID(), indx);  break;
  case BC_CFL:  this->connectivity_assembling(false);  ((Geometry*)Geom)->assign_FaceLoad_by_elem_ID(bc, this->give_ID(), indx);  break;
  case BC_CBL:  
  case BC_DW:                                          ((Geometry*)Geom)->assign_BodyLoad_by_elem_ID(bc, this->give_ID()      );  break;
  default: _errorr("unsupported bc type");
  }
}


void Element :: set_mprop (long val)
{
  GeometryComponent :: set_mprop(val);
  elemAttribs()->set_Patt();
}
void Element :: set_prop_node_inher (bool everynode)
{
  long i;

  if (this->give_mproperty_cnt() == 0)  errol;

  for (i=0; i<points(); i++)
    if (everynode || points[i]->give_numsuperelem() == 1)
      ((Point*)points[i])->set_mprop(this->give_mproperty());

}

void Element :: integrate_duplicated_cell (const Element *slave)
{
  Cell :: integrate_duplicated_cell (slave);

  elemAttribs()->integrate_duplicated_one(slave->give_elemAttribs());
}




void Element :: switch_node_pointer_in_all_components (Point *slave, Point *master, bool duplcheck)
{
  this->switch_node_pointer(slave, master, duplcheck);
  long i;
  for (i=0; i<edges(); i++)  ((Edge*)edges[i])->switch_node_pointer(slave, master, duplcheck);
  for (i=0; i<faces(); i++)  ((Face*)faces[i])->switch_node_pointer(slave, master, duplcheck);
}

void Element :: switch_node_pointer (Point *slave, Point *master, bool duplcheck)
{
  Cell :: switch_node_pointer (slave, master, duplcheck);

  if (duplmaster == this  ||  !duplcheck)  return;

  for (int i=0; i<master->give_numsuperelem(); i++)
    if ( master->give_superelem(i)->has_same_geom_with(this) ) {
      if ( this->give_ID() > master->give_superelem(i)->give_ID() )
    this->setup_duplicity_master((Cell*)master->give_superelem(i));
      else
    ((Cell *)master->give_superelem(i))->setup_duplicity_master(this);
    }
}


bool Element :: is_point_on (const PoinT *point, const GeometryComponent *&comp, PoinT *nc) const
{
  // check, jestli je point vubec v circum
  if ( ! const_cast<Element*>(this)->is_point_in_sphere(point) )
    return false;

  long i;

  // check, jestli se point nerovna nejakemu nodu
  for (i=0; i<points(); i++)
    if ( points[i]->is_identical_to (this->give_circum(), point) ) {
      comp = points[i];
      return true;
    }

  // check, jestli neni point nahodou on an edge
  for (i=0; i<edges(); i++)
    if (edges[i]->is_nod_on (this->give_circum(), point, nc->x)) {
      comp = edges[i];
      return true;
    }

  // check, jestli neni point nahodou na ??
  const Face *fac=NULL;
  const Edge *t1, *t2;  t1 = t2 = NULL;

  switch (this->give_classid()) {
  case classBrick:
    double x[8],y[8],z[8];
    for (i=0;i<8;i++) {
      x[i] = this->give_point(i)->give_coord(0);  //                                                           *****
      y[i] = this->give_point(i)->give_coord(1);  //                                                           *(2)*
      z[i] = this->give_point(i)->give_coord(2);  //                                                           *****
    }
    if (nc_brick_3d (ZERO,x,y,z,point,nc)) {
      if      ( fabs(nc->x-1.0) < ZERO ) { fac=faces[4]; nc->x=nc->y; nc->y=nc->z; t1=edges[10]; t2=edges[ 2]; } // ksi =  1.0
      else if ( fabs(nc->x+1.0) < ZERO ) { fac=faces[2]; nc->x=nc->y; nc->y=nc->z; t1=edges[ 9]; t2=edges[ 0]; } // ksi = -1.0    *****
      else if ( fabs(nc->y-1.0) < ZERO ) { fac=faces[3];              nc->y=nc->z; t1=edges[10]; t2=edges[ 1]; } // eta =  1.0    *(2)*
      else if ( fabs(nc->y+1.0) < ZERO ) { fac=faces[5];              nc->y=nc->z; t1=edges[11]; t2=edges[ 3]; } // eta = -1.0    *****
      else if ( fabs(nc->z-1.0) < ZERO ) { fac=faces[0];                           t1=edges[ 2]; t2=edges[ 1]; } // dzeta =  1.0
      else if ( fabs(nc->z+1.0) < ZERO ) { fac=faces[1];                           t1=edges[ 6]; t2=edges[ 5]; } // dzeta = -1.0
      else {
    // volume
    comp = this;
    return true;
      }

      // face
      fac->transform_nc (t1, t2, nc);
      comp = fac;
      return true;

    }
    else return false;

    break;
  default: _errorr("Unknown element type is required in function is_point_on_el");
  }

  return false;
}

void Element :: print_row_VTK (FILE *stream) const
{
  // nodes
  fprintf (stream, "%ld", points());
  for (int i=0; i<points(); i++)
    fprintf (stream, " %ld", ((Node*)points[i])->give_lid_id(0));

  fprintf (stream, "\n");
}
void Element :: print_row_VTX (char *str) const
{
  // nodes
  str[0]='\0';
  for (int i=0; i<points(); i++)
    sprintf (strchr(str, '\0'), " %ld", ((Node*)points[i])->give_lid_id(0));
}



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS GELEMENT                        *** *** *** ***
//*  ********************************************************************************************
Gelement :: Gelement (classID mecg, long gid, long oid, const Geometry *mg, long nn, long ne, long nf) : Element (mecg, gid, oid, 0, mg, nn, ne, nf)
{
  attributes = new GelemAttribs(this, (long)0, EAL_direct);
  fixedVertices = NULL;
  fixedGelements = NULL;
}
Gelement :: Gelement (classID mecg, const Gelement * src) : Element(mecg, src)
{
  //
  attributes = new GelemAttribs(this, (const GelemAttribs*)src->attributes, EAL_direct);
}

void Gelement :: set_elemSize (double val)
{
  if (this->give_nfa()) {
    if (this->give_nfa() != 1)  _errorr("elemSize at element with more than 1 face");
    
    ((Face*)this->give_face(0))->give_facedgeAttribs()->set_elemSize(val);
  }
  else {
    if (this->give_ned() != 1)  _errorr("elemSize at element with more than 1 edge");
    
    ((Edge*)this->give_edge(0))->give_facedgeAttribs()->set_elemSize(val);
  }
}

void Gelement :: set_elemCount (double val)
{
  if (this->give_ned() != 1)  _errorr("elemCount at element with more than 1 edge");
  
  ((Edge*)this->give_edge(0))->give_facedgeAttribs()->set_elemCount(val);
}

void Gelement :: assign_fixed_entities_by_ID (bool node, long ncn, const long *icn)
{
  if (node) { if (fixedVertices  == NULL)  fixedVertices  = new GPA<const Vertex>;   }
  else      { if (fixedGelements == NULL)  fixedGelements = new GPA<const Gelement>; }
  
  if (node) for (long i=0; i<ncn; i++)  fixedVertices ->add(((Model*)Geom)->give_Vertex  (icn[i]-1));
  else      for (long i=0; i<ncn; i++)  fixedGelements->add(((Model*)Geom)->give_Gelement(icn[i]-1));
}


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS POLYLINE                        *** *** *** ***
//*  ********************************************************************************************
void PolyLine :: checkConsistency (void) const
{
  Gelement :: checkConsistency();
  
  //* Detects (poly)line with hinges/condensed rotation DOFs at both end point.
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_Dxy_R___:
  case DPN_DxyzR___: break;
  case DPN_Dxy_R__z:
  case DPN_DxyzRxyz:
    if (this->elemAttribs()->has_fullhinge_ends())
      _errorr2("Rotation DOFs condensed at both ends of Polyline - forbidden kinematic mechanism at elemnet %ld.", ID+1);
    break;
  default: errol;
  }
}

void PolyLine :: init_point_on (Mesh* pm, GPA<const Node> &hnodes, const PoinT *point, const GeometryComponent* comp, const PoinT *nc, const Point *parentpoint) const
{
  const Node* node = dynamic_cast<const Node *> (comp);
  if (node) {
    hnodes.add(node);
    return;
  }

  const Cell* cell = dynamic_cast<const Cell *> (comp);  if (!cell) errol;

  HangingNode *hn = new HangingNode(pm, (long)-1);

  int ord = -4;
  //
  const Facedge *faed = dynamic_cast<const Facedge*>(cell);
  if (faed)
    ord = faed->give_ord();
  else {
    const FElement *elem = dynamic_cast<const FElement*>(cell);
    if (elem)
      ord = elem->give_ord();
    else errol;
  }

  // !!!!!!!!!!!!!!
  hn->initialize_hn (point, cell, cell->give_dimension(), ord, cell->give_nno(), (const Node**)cell->give_points()->v(), nc);
  // set property
  if (parentpoint) { if (parentpoint->give_mproperty_cnt())  hn->set_mprop(parentpoint->give_mproperty()); }
  // tady hn zdedil property z elementu, byla prirazena do mproperty, to jsem ale zmenil:
  //else             { if (       this->give_mproperty_cnt())  hn->set_mprop(       this->give_mproperty()); }
  // zmenil jsem to na prirazeni do vproperty
  if (       this->give_mproperty_cnt())  hn->add_property (3, this->give_mproperty());
  
  hnodes.add((const Node*)hn);
  ((Mesh*)pm)->add_another_Pjnt ((Point*)hn);
}

PolyLine* PolyLine :: generate_mesh_RFbyHN(Mesh* pm, GPA<const Node> &hnodes) const
{
  long i;
  PolyLine *pol1, *pol2 = NULL;

  if (hnodes()==0) {
    const GeometryComponent *comp;
    PoinT nc;  nc.zero();
    int stpid = 0;
    while (true) {
      if ( pm->give_3Delement_with_point_inside (points[stpid]->give_coords(), comp, &nc) ) {
    // first point is inside
    if (stpid == 0)
      this->init_point_on (pm, hnodes, points[0]->give_coords(), comp, &nc, points[0]);
    // last point is inside => reverse polyline
    else if (stpid == points()-1) {
      ((PolyLine*)this)->points.reverse();
      this->init_point_on (pm, hnodes, points[0]->give_coords(), comp, &nc, points[0]);
    }
    // middle point is inside => split polyline
    else {
      pol1 = (PolyLine*)this;
      pol2 = new PolyLine(this);

      pol2->points.reverse();
      pol2->del_points_down(points()-stpid);
      pol2->points.reverse();

      pol1->del_points_down(1       +stpid);
      pol1->points.reverse();

      this->init_point_on (pm, hnodes, points[0]->give_coords(), comp, &nc, points[0]);
    }

    break;
      }
      else {
    if (Pd->give_PDBO(PDBO_cutRF)) {
      // first point is outside => try the last
      if      (stpid == 0)           stpid = points()-1;
      // last point is outside => try the second
      else if (stpid == points()-1)  stpid = 1;
      // all point are outside => warning
      else if (stpid == points()-2) { _warningg2("Entire polygon (id %ld) lies outside of primary mesh", this->give_ID()); }
      // n-th point is outside => try n+1
      else     stpid++;
    }
    else
      // first point is outside => error
      _errorr2 ("Start point of RF bar %ld doesn't lay inside of domain", this->give_ID());
      }
    }
  }

  this->divide(pm, hnodes);

  if (pm->give_Parallel()) {
    HangingNode *hn;

    for (i=0; i<hnodes(); i++) {
      hn = dynamic_cast<HangingNode *>((Node*)hnodes[i]);
      if (hn)  ((HNAttribs*)hn->give_pointAttribs())->find_hndomain ();
    }
  }

  this->findout_segment_domain (pm, hnodes);

  return pol2;
}

void PolyLine :: divide (Mesh* pm, GPA<const Node> &hnodes) const
{
  int dimMC;
  long i,j,tc,konec,clastelem;
  long csuperel=0,cunwnod=0,cunwfac=0; //
  PoinT ip, nc;
  const PoinT *curvert;
  const Node *nod, *curnod;
  const Point *auxnod, **unwnod=NULL;
  const HangingNode *curhnod;
  const Cell *compid;
  const Face    *face, **unwfac=NULL;                  // unwanted faces
  const Edge    *edge;
  const Element *elem, **superel=NULL, **lastelem;
  const GeometryComponent *comp;
  
  nc.zero();
  
  // loop over all segments
  for (i=1; i<points(); i++) {
    curvert = points[i]->give_coords();
    clastelem = 0;
    lastelem = NULL;
    konec = 0;
    
    // it goes along polyline to the end (konec) point
    while (!konec) {
      curnod  = hnodes.last();
      curhnod = dynamic_cast<const HangingNode *> (curnod);
      
      if (curhnod) dimMC = curhnod->give_HNattrb()->give_dimMC();
      else         dimMC = 0;
      
      // najde superior elements, kde by teoreticky mohl; byt dalsi HN
      // find superior elements = elements with potential HN
      switch (dimMC){
      case 0:
    nod = curnod;
    
    csuperel  = nod->give_numsuperelem();          superel = nod->give_superelems()->v();
    cunwnod = 1;                                   unwnod  = (const Point**)&nod;
    cunwfac = nod->give_numsuperface();            unwfac  = nod->give_superfaces()->v();
    break;
      case 1:
    edge = (const Edge*)curhnod->give_HNattrb()->give_MC();
    
    csuperel = edge->give_numsuperelem();          superel = edge->give_superelems()->v();
    cunwfac  = edge->give_numsuperface();          unwfac  = edge->give_superfaces()->v();
    if (clastelem) { cunwnod = edge->give_nno();  unwnod = edge->give_points()->v(); }
    else           { cunwnod = 0;                 unwnod = NULL;       }
    break;
      case 2:
    face = (const Face*)curhnod->give_HNattrb()->give_MC();
    
    csuperel = face->give_numsuperelem();          superel = face->give_superelems()->v();
    cunwfac = 1;                                   unwfac  = &face;
    if (clastelem) { cunwnod = face->give_nno(); unwnod = face->give_points()->v(); }
    else           { cunwnod = 0;                unwnod = NULL;       }
    break;
      case 3:
    elem = (const Element*)curhnod->give_HNattrb()->give_MC();
    
    csuperel = 1;                                  superel  = (const Element**)&elem;
    cunwnod = 0;                                  unwnod = NULL;
    cunwfac = 0;                                  unwfac = NULL;
    break;
      default:  _errorr ("Unknown dimension of master node");
      }

      // find intersection point of "superior element" and polyline
      for (j=0; j<csuperel; j++){
    if (same_array_elements_asym (1,superel+j,clastelem,lastelem))   continue;
    if (superel[j]->give_dimension() != 3)  continue;

    // segment konci inside of element
    if (((FElement*)superel[j])->is_point_on (curvert, comp, &nc)) {
      this->init_point_on (pm, hnodes, curvert, comp, &nc, points[i]);

      konec = 1;
      break;
    }

    // segment protina node of element
    if (((FElement*)superel[j])->cross_abscissa_node (curnod->give_coords(), curvert, cunwnod, unwnod, auxnod, nc.x, &ip)){
      this->init_point_on (pm, hnodes, NULL, auxnod, &nc, NULL);
      break;
    }

    // segment intersects edge or face of element
    tc = ((FElement*)superel[j])->cross_abscissa_face (curnod->give_coords(), curvert, cunwfac, unwfac, compid, &nc, &ip);

    if (tc==1 || tc==2) { this->init_point_on (pm, hnodes, &ip, compid, &nc, NULL);  break; }
      }

      if (j==csuperel) {
    if (Pd->give_PDBO(PDBO_cutRF)) {
      if (dimMC==2  &&  csuperel==1  &&  superel[0]->give_dimension()==3) {
        konec = 1;
        i = points();
      }
      else {
        if (dimMC==2  &&  csuperel==2  &&  superel[0]->give_dimension()==3)
          errol;
        else
          errol;
      }
    }
    else
      _errorr ("Something is wrong in function polyline_dividing");
      }

      lastelem = superel;
      clastelem = csuperel;

    } // end of while

  } // end of for

  //for (i=0;i<Np;i++){
  //  if (PolyLines[i]->te==LiBeam){
  //    for (j=0;j<PolyLines[i]->nhn;j++){
  //    if (PolyLines[i]->hnodes[j]->dimMC==0)
  //      Nodes[ PolyLines[i]->hnodes[j]->idMC ]->property = 7;
  //    else
  //      fprintf (stderr,"\n\n unknown type of dimmension of master component is required in function findout_hnode_domain. (%s, line %d).\n",__FILE__,__LINE__);
  //    }
  //  }
  //}
}

void PolyLine :: findout_segment_domain (Mesh* pm, GPA<const Node> &hnodes) const
{
  long i,j,dom=0;
  FElement *elem = NULL;

  for (i=1; i<hnodes(); i++) {
    // find domain
    if (!pm->give_Parallel()) dom = 0;
    else
      if      (!hnodes[i-1]->is_shared())  dom = hnodes[i-1]->give_domain();
      else if (!hnodes[i  ]->is_shared())  dom = hnodes[i  ]->give_domain();
      else
    for (j=0; j<pm->give_NumDomains(); j++)
      if (hnodes[i-1]->give_lid(j)>=0 && hnodes[i  ]->give_lid(j)>=0)
        { dom = j; break; }
    
    // allocate new segment
    if (Pd->give_fulldata())
      if (this->elemAttribs()->give_dpn() != DPN_DxyzR___)  _errorr ("Unknown type of segment");
    
    elem = new Beam (-1, -1, pm, this, dom, -1);
    
    elem->set_node(0, hnodes[i-1]);
    elem->set_node(1, hnodes[i  ]);
    elem->set_mprop(this->give_mproperty());

    pm->add_another_Element (elem, dom);
  }
}


//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS LINE                          *** *** *** ***
//*  ********************************************************************************************
void Line :: initialize (void)
{
  PolyLine :: initialize ();
  
  //* The line with hinges/condensed rotation DOFs at both end point is forced to be truss type.
  //* Tady by mela byt kontrola, ze na obou koncich neni uvolnena rotace kolem x,
  //* pokud je, tak predepisu zakaz rotace kolem osy x, ale jen prvnimu uzlu!!
  //* prerod na truss muzu udelat jen pokud argument prikaz radky, a pokud neni stabilita a pokud...
  //* Az to naimplementujes, tak rozchod udrzovane/stability/data/TEST_dodelat
  if (Pd->give_PDBO (PDBO_hingedBeamToTruss)) {
    if (this->give_nno() != 2)  errol;
    
    this->elemAttribs()->switch_dpn_Line();
  }
}
long Line :: give_edge_nodes (const Point **&edgnodes) const
{
  if (points() != 2) errol;
  edgnodes = new const Point*[2];

  edgnodes[0] = points[0];
  edgnodes[1] = points[1];
  return 2;
}


//*  ********************************************************************************************
//*  *** *** *** ***                       CLASS POLYGONMDL                       *** *** *** ***
//*  ********************************************************************************************
void PolygonMdl :: allocNnNe (int n)
{
  ((ComponentGeometry2Dpolygon*)cg)->set_nv(n);
  points.resizeup(n);
  edges. resizeup(n);
}

long PolygonMdl :: give_edge_nodes (const Point **&edgnodes) const
{
  edgnodes = new const Point*[2*edges()];

  int i;
  for (i=0; i<edges()-1; i++) {
    edgnodes[2*i  ] = points[i  ];
    edgnodes[2*i+1] = points[i+1];
  }
  edgnodes[2*i  ] = points[i];
  edgnodes[2*i+1] = points[0];
  return 2;
}
long PolygonMdl :: give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const
{
  facnodes = new const Point* [points()];
  facedges = new const Edge*  [points()];

  for (int i=0; i<points(); i++) {
    facnodes[i] = points[i];
    facedges[i] = edges[i];
  }
  return points();
}



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS FELEMENT                        *** *** *** ***
//*  ********************************************************************************************
FElement :: FElement (classID mecg, long gid, long oid, const Geometry *mg, long o, long nn, long ne, long nf, bool aa, long dom, long lid) : Element (mecg, gid, oid, 0, mg, nn, ne, nf)
{
  if (aa)  this->attributes_allocation(NULL);
  
  //
  FElement::domain = dom;
  FElement::lid = lid;
  
  mdl_masterel = NULL;
  
  //
  results = NULL;
  maxSigmaEq = NULL;
  regid = -1;
}
FElement :: FElement (classID mecg, const FElement * src) : Element(mecg, src)
{
  //
  attributes = new FElemAttribs(this, (const FElemAttribs*)src->attributes, '?');

  domain = 0;        if (src->domain)          errol;
  lid = ID;          if (src->lid != src->ID)  errol;
  results = NULL;    if (src->results)         errol;
  maxSigmaEq = NULL; if (src->maxSigmaEq)      errol;
}
FElement :: ~FElement()
{
  deallocateCheckUno  (results, Msh()->give_rslts_nsteps(), cRTE, false);
  if (maxSigmaEq)  delete maxSigmaEq;
}

void FElement :: attributes_allocation (const GelemAttribs *masterat)
{
  if (attributes)  errol;
  
  if (masterat)  attributes = new FElemAttribs(this, masterat, 'm');
  else           attributes = new FElemAttribs(this, (long)0, EAL_direct);
}


void FElement :: initialize (void)
{
  Element :: initialize ();
}
//void FElement :: finitialize (void){  Element :: finitialize ();  //attributes->finitialize();}
void FElement :: checkConsistency (void) const
{
  Element :: checkConsistency();
  
  if (elemAttribs()->give_ord() > 0)
    _errorr("odkomentuj nasledujici");
  // bool ordcheck;
  // if ( Msh()->connectivity_is_assembled() )  ordcheck = true;
  // else                                          ordcheck = false;
  // 
  // if (ordcheck) {
  //   long i;
  // 
  //   // *** check same ord ***
  //   for (i=0; i<edges(); i++)
  //     if (elemAttribs()->give_ord() != edges[i]->give_ord())
  //    _errorr2 ("ord != edges[%d]->ord", i);
  // 
  //   // *** check same ord ***
  //   for (i=0; i<faces(); i++)
  //     if (elemAttribs()->give_ord() != faces[i]->give_ord())
  //    _errorr2 ("ord != faces[%d]->ord", i);
  // }
}

void FElement :: set_model_prop (long val, const Model *model, bool flag)
{
  GeometryComponent :: set_mprop (val);
  
  mdl_masterel = dynamic_cast<const Gelement*>(model->give_Elem(val-1));
  
# ifdef DEBUG
  if (mdl_masterel == NULL)  errol;
  if (this->give_dimension() !=  mdl_masterel->give_dimension())  errol;
  if (flag)  errol;
# endif
}
 

IntPointSet FElement :: give_IPset_comp (Solver sol) const
{
  IntPointSet IPset = this->elemAttribs()->give_IPset();
  if (IPset != IPS_Void)  return IPset;
  
  FiniteElementTypeSet fetset;
  return IntPointSet_fet2e_comp (FETSet_set2e (this->elemAttribs()->give_FETS(&fetset), sol, Pd->give_analgroup()));
  
  return IPS_Void;
}

IntPointSet FElement :: give_IPset_rslts (Solver sol) const
{
  FiniteElementTypeSet fetset;
  return IntPointSet_fet2e_rslts (this->give_IPset_comp(sol), FETSet_set2e (this->elemAttribs()->give_FETS(&fetset), sol, Pd->give_analgroup()));
}

// /// higher set = for displacement^2 computation, or mass matrix or sigma error ...
// IntPointSet FElement :: give_IPset2 (void) const
// {
//   //IntPointSet IPset; // = this->elemAttribs()->give_IPset2();
//   
//   /// default IPset
//   FiniteElementTypeSet fetset;
//   return IntPointSet_fets2ips_2 (this->elemAttribs()->give_FETS(&fetset));
// }


void FElement :: integrate_duplicated_cell (const Element *slave)
{
  Element :: integrate_duplicated_cell (slave);
  
  FiniteElementTypeSet fets1;  this->            elemAttribs()->give_FETS(&fets1);
  FiniteElementTypeSet fets2;  ((FElement*)slave)->elemAttribs()->give_FETS(&fets2);
  
  if (! fets1.is_equal_to(fets2))  _errorr ("Duplicity slave and master differ in rot atribut");
  
  const FElement *slv = dynamic_cast<const FElement *>(slave);
  if (!slv)  _errorr ("The slave is not same type as this");
}

 // const Dvctr* FElement :: give_results_d (long step, ResultTypesAtElem rt) const
 // {
 //   if (results[step][rt] == NULL)  return NULL;
 //
 //   Dvctr *rslt = dynamic_cast<Dvctr*>(results[step][rt]);
 //   if (!rslt) errol;
 //
 //   return rslt;
 // }


//*  *******************
//*  ***   RESULTS   ***
//*  *******************
void FElement :: allocate_results (void)
{
  if (results) return;

  results = new Array**[Msh()->give_rslts_nsteps()];

  for (long i=0; i<Msh()->give_rslts_nsteps(); i++) {
    results[i] = new Array*[cRTE];
    memset (results[i], 0, cRTE*sizeof(Array*));
  }
}


void FElement :: add_result (         Array *rslt, long step, ResultTypesAtElem rt) { this->allocate_results();  if (this->results[step][rt])  errol;  this->results[step][rt] = rslt;}
void FElement :: set_result (long s, double *rslt, long step, ResultTypesAtElem rt) { this->allocate_results();  if (this->results[step][rt])  errol;  this->results[step][rt] = new Dvctr(s, rslt); }
void FElement :: set_result (        double  rslt, long step, ResultTypesAtElem rt) { this->allocate_results();  if (this->results[step][rt])  errol;  this->results[step][rt] = new Dscal(   rslt); }
void FElement :: set_result (  const VectoR *rslt, long step, ResultTypesAtElem rt) { this->allocate_results();  if (this->results[step][rt])  errol;  this->results[step][rt] = new Dvctr(   rslt); }

const Dscal* FElement :: give_results_ds (long step, ResultTypesAtElem rt) const { if (results[step][rt] == NULL)  return NULL;  Dscal *rslt = dynamic_cast<Dscal*>(results[step][rt]);  if (!rslt) errol;  return rslt; }
const Dvctr* FElement :: give_results_dv (long step, ResultTypesAtElem rt) const { if (results[step][rt] == NULL)  return NULL;  Dvctr *rslt = dynamic_cast<Dvctr*>(results[step][rt]);  if (!rslt) errol;  return rslt; }
const Dmtrx* FElement :: give_results_dm (long step, ResultTypesAtElem rt) const { if (results[step][rt] == NULL)  return NULL;  Dmtrx *rslt = dynamic_cast<Dmtrx*>(results[step][rt]);  if (!rslt) errol;  return rslt; }

void FElement :: setup_maxSigmaEq (void)
{
  if (maxSigmaEq)  return;
  
  int nl = this->elemAttribs()->give_cs()->give_nlayers();
  
  if (nl) maxSigmaEq = new double[nl];
  else    maxSigmaEq = new double;
  
  this->compute_maxSigmaEq();
  
  if (maxSigmaEq[0] <= 0.0) {
    delete maxSigmaEq;
    maxSigmaEq = NULL;
  }
}

double FElement :: give_CSusage_elast ()
{
  this->setup_maxSigmaEq();
  
  if (maxSigmaEq == NULL)  return 0.0;
  
  int nl = this->elemAttribs()->give_cs()->give_nlayers();
  
  if (nl == 0)  return maxSigmaEq[0];
  
  double max = 0.0;
  
  for (int i=0; i<nl; i++)
    if (max < maxSigmaEq[i])
      max = maxSigmaEq[i];
  
  return max;
}

double FElement :: give_CSusage_elast_rel ()
{
  this->setup_maxSigmaEq();
  
  if (maxSigmaEq == NULL)  return 0.0;
  
  int nl = this->elemAttribs()->give_cs()->give_nlayers();
  
  if (nl == 0)  return maxSigmaEq[0] / elemAttribs()->give_mat()->give_Ry();
  
  double val, max = 0.0;
  
  for (int i=0; i<nl; i++) {
    val = maxSigmaEq[i] / this->elemAttribs()->give_cs()->give_mat(i)->give_Ry();
    if (max < val)
      max = val;
  }
  
  return max;
}
int FElement :: give_CSusage_elast_bool ()
{
  double answer = this->give_CSusage_elast_rel();
  
  return ( answer > 0.0 ? answer > 1.0 : -1 );
}

double FElement :: fillupbyzero (Dvctr *data, SStype SST) const
{
  switch (SST) {
  case SST_truss:     data->resize_ignore_vals(1);                            data->zero();  break;
  case SST_plstrain:
  case SST_plstress:  data->resize_ignore_vals(3);                            data->zero();  break;
  case SST_beam:      data->resize_ignore_vals(9);                            data->zero();  break;
  case SST_3d:
  case SST_3dshell:   data->resize_ignore_vals(6);                            data->zero();  break;
  case SST_transp2d:  data->resize_ignore_vals(3*Pd->give_global_nDOFs());  data->zero();  break;

  default:                _errorr2("unsupported SST_* %d", SST);
  }

  return 0.0;
}

long FElement :: give_result_ncomp (long time_step, ResultTypesAtElem rte) const
{
  if (results == NULL)  errol;
  
  if (results[time_step][rte] == NULL) return 0;
  
  Xmtrx *mtrx = dynamic_cast<Xmtrx*>(results[time_step][rte]);
  if (mtrx == NULL)  errol;
  
  return mtrx->give_ccols();
}

void FElement :: read_input (const char *&str, femFileFormat fff)
{
  switch (fff) {
  case FFF_T3d: {
    long parentdim, auxl, prop;
    
    this->read_nodes (str, fff);
    
    SP_scan_number (str, parentdim);
    SP_scan_number (str, auxl);
    SP_scan_number (str, prop);
    
    // master model element
    const Model *model = ((Mesh*)Geom)->give_Master_model();
    
# ifdef DEBUG
    if (model && prop<1)  errol;
# endif
    
    if (model) {
      this->set_model_prop (prop, model);
      this->attributes_allocation(mdl_masterel->give_gelemAttribs());
    }
    else {
      this->attributes_allocation(NULL);
      this->set_mprop(prop); // attributes allocation have to be before property setting
    }
    
    this->connectivity_assembling(false);
    
    switch (this->give_cellGeom()) {
    case CG_Line:       if (model)  ((Edge*)edges[0])->set_model_prop (prop, model, true);  break;
    case CG_Triangle:
    case CG_Quadrangle: if (model)  ((Face*)faces[0])->set_model_prop (prop, model, true);
      // edge property
      SP_skip_word (str, edges());
      for (int i, e=0; e<edges(); e++) {
    i = ECN_convert_e2i (e, 1, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
    SP_scan_number (str, prop);
    if (prop) {
      if (model) ((Edge*)edges[i])->set_model_prop (prop, model, false);
      else       ((Edge*)edges[i])->checkset_mprop(prop);
    }
# ifdef DEBUG
    else if (    ((Edge*)edges[i])->give_mproperty_cnt())  errol;
# endif
      }
      break;
    case CG_Hexahedron:
      // face property
      for (int i, e=0; e<faces(); e++) {
    i = ECN_convert_e2i (e, 2, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
    SP_scan_number (str, prop);
    if (prop) {
      if (model) ((Face*)faces[i])->set_model_prop (prop, model, false);
      else       ((Face*)faces[i])->checkset_mprop(prop);
    }
# ifdef DEBUG
    else if (    ((Face*)faces[i])->give_mproperty_cnt())  errol;
# endif
      }
      break;
    default:  _errorr2("unsupported cell geometry CG_* %d", this->give_cellGeom());
    }
    
    break;
  }
  case FFF_OOFEM:
    SP_scan_expected_word_exit   (str, "nodes", "Invalid structure of the given OOFEM input file", CASE);
    SP_scan_expected_number_exit (str, points(), "Invalid input, number of nodes on element");
    
    this->read_nodes (str, fff);
    
    //if (attributes)  _errorr ("attributes already allocated");
    //attributes = new ElemAttribs(0, Pd, this);
    elemAttribs()->initialize_from (str, fff);
    break;
  case FFF_SIFEL:
    this->read_nodes (str, fff);
    
    int auxi;
    SP_scan_number (str, auxi);  if (auxi) errol; // code numbers
    
    elemAttribs()->initialize_from (str, FFF_SIFEL);
    break;
  case FFF_JKTK:
    this->read_nodes (str, fff);
    
    this->connectivity_assembling(false);
    
    long prop;
    SP_scan_number (str, prop);
    this->set_mprop(prop);
    
    int e, i;
    for (e=0; e<edges(); e++) {
      SP_scan_number (str, prop);
      i = ECN_convert_e2i (e, 1, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      if (prop) ((Edge*)edges[i])->checkset_mprop(prop);
# ifdef DEBUG
      else if ( ((Edge*)edges[i])->give_mproperty_cnt())  errol;
# endif
    }
    for (e=0; e<faces(); e++) {
      SP_scan_number (str, prop);
      i = ECN_convert_e2i (e, 2, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      if (prop) ((Face*)faces[i])->checkset_mprop(prop);
# ifdef DEBUG
      else if ( ((Face*)faces[i])->give_mproperty_cnt())  errol;
# endif
    }
    break;
  default: _errorr("unsupported fem file format");
  }
}


void FElement :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  int e;
  long i;
  FiniteElementTypeSet fets;
  
  switch (fff) {
  case FFF_OOFEM:
    // name
    fprintf (stream, "%s", FETSet_set2s_OOFEM (this->elemAttribs()->give_FETS(&fets), Pd->give_analgroup()));
    
    // nodes
    if      (Parallel())                           fprintf (stream, " %4ld  nodes %ld  ",  lid+1, points());
    else if (Pd->give_PDBO (PDBO_OUT_origelemid))  fprintf (stream, " %4ld  nodes %ld  ", origid, points());
    else                                           fprintf (stream, " %4ld  nodes %ld  ",   ID+1, points());
    
    for (e=0; e<points(); e++) {
      i = ECN_convert_e2i (e, 0, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      fprintf (stream, "%4ld ", ((Node*)points[i])->give_lid_id(domain)+1);
    }
    
    // attributes
    elemAttribs()->print_row (stream, fff, domain);
    break;
  case FFF_SIFEL:
    // nodes
    fprintf (stream, "%6ld", this->ID+1);
    fprintf (stream, "  %d", FETSet_set2i_SIFEL (this->elemAttribs()->give_FETS(&fets), Pd->give_analgroup()));
    if (elemAttribs()->give_sst() == SST_plstress)  fprintf (stream, "  10");
    if (elemAttribs()->give_sst() == SST_plstrain)  fprintf (stream, "  11");
    
    for (e=0; e<points(); e++) {
      i = ECN_convert_e2i (e, 0, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      fprintf (stream, " %6ld", points[i]->give_ID()+1);
    }
    
    // attributes
    elemAttribs()->print_row (stream, fff, domain);
    break;
  case FFF_ANSYS:
    // type
    fprintf (stream, "TYPE, %d\n", FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_ANSYS, Pd->give_analgroup()) - FET_ANSYS_start);
    
    // attributes
    elemAttribs()->print_row (stream, fff, domain);
    
    fprintf (stream, "E");
    for (e=0; e<points(); e++) {
      i = ECN_convert_e2i (e, 0, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      fprintf (stream, ",%6ld", points[i]->give_ID()+1);
    }
    
    fprintf (stream,"\n");
    break;
  case FFF_JKTK: {
    bool m = Pd->give_PDBO(PDBO_melnik);
    
    // nodes
    if (!m) fprintf (stream, "%8ld", this->ID+1);
    else    fprintf (stream, "%ld", this->ID+1);
    fprintf (stream, " %d", CellGeometry_e2i_JKTK (this->give_cellGeom()));
    if (!m) fprintf (stream, "  ");
    
    for (e=0; e<points(); e++) {
      i = ECN_convert_e2i (e, 0, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      if (!m) fprintf (stream, "%7ld", points[i]->give_ID()+1);
      else    fprintf (stream, " %ld", points[i]->give_ID()+1);
    }
    
    if (!m) fprintf (stream, "%5ld", this->give_mpropertyORzero());
    else    fprintf (stream, " %ld", this->give_mpropertyORzero());
    
    if (!m) { if (edges())  fprintf (stream, " "); }
    for (e=0; e<edges(); e++) {
      i = ECN_convert_e2i (e, 1, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      fprintf (stream, " %ld", edges[i] ? edges[i]->give_mpropertyORzero() : 0);
    }
    if (!m) { if (faces())  fprintf (stream, " "); }
    for (e=0; e<faces(); e++) {
      i = ECN_convert_e2i (e, 2, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
      fprintf (stream, " %ld", faces[i] ? faces[i]->give_mpropertyORzero() : 0);
    }
    break;
  }
  default: errol;
  }
  
  //
  if (endline)  fprintf (stream,"\n");
}

void FElement :: read_nodes (const char *&str, femFileFormat fff)
{
  int e, i;
  long auxl;
  for (e=0; e<points(); e++) {
    SP_scan_number (str, auxl);
    i = ECN_convert_e2i (e, 0, this->elemAttribs()->give_ord(), this->give_cellGeom(), fff);
    this->set_node(i, auxl-1);
  }
}

bool FElement :: is_cmfr (void)
{
//  docasne zakomentovano  if (nno != 2) _errorr ("nno != 2") ;
//  docasne zakomentovano
//  docasne zakomentovano  const Node *node1 = nodes[0];   if ( node1->give_pointAttribs()->give_BC()[0] ) return false;
//  docasne zakomentovano  const Node *node2 = nodes[1];   if ( node2->give_pointAttribs()->give_BC()[0] ) return false;
//  docasne zakomentovano
//  docasne zakomentovano  long nexn1 = node1->give_numsuperelem() - 1;    //if (nexn1 == 0)  return true;
//  docasne zakomentovano  long nexn2 = node2->give_numsuperelem() - 1;    //if (nexn2 == 0)  return true;
//  docasne zakomentovano
//  docasne zakomentovano  // smaze pruty ktery se dotykaji nepodepreneho uzlu ktereho se dotykaji jen 1 nebo 2 pruty
//  docasne zakomentovano  if ( nexn1<2 || nexn2<2 ) {
//  docasne zakomentovano    this->set_delete_flag(true);
//  docasne zakomentovano    return true;
//  docasne zakomentovano  }

  return false;

  // smaze prouty ktery netvori trojuhelnik s jinymi
  // u vlny je to bohuzel 100%
  //
  //const Node **exns = new const Node* [nexn1+nexn2];   nexn1 = 0;
  //
  //for (int i=0; i<node1->give_numsuperelem(); i++)
  //  if (node1->give_superelem(i) != this) {
  //    if (node1->give_superelem(i)->give_node(0) != node1)  exns[nexn1++] = node1->give_superelem(i)->give_node(0);
  //    else                                                exns[nexn1++] = node1->give_superelem(i)->give_node(1);
  //  }
  //
  //for (int i=0; i<node2->give_numsuperelem(); i++)
  //  if (node2->give_superelem(i) != this) {
  //    if (node2->give_superelem(i)->give_node(0) != node2)  exns[nexn1++] = node2->give_superelem(i)->give_node(0);
  //    else                                                exns[nexn1++] = node2->give_superelem(i)->give_node(1);
  //  }
  //
  //nexn2 = members_are_unique (nexn1, exns);
  //
  //delete [] exns;
  //
  //if (nexn2) {
  //  this->set_delete_flag(true);
  //  return true;
  //}
  //
  //return false;
}


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS ELEMENT                         *** *** *** ***
//*  ********************************************************************************************
void OOFEM_output_scan_elem_head (const char *str, long o, long i, bool oid)
{
  ;        SP_scan_expected_word_exit   (str, "element", "Invalid structure of OOFEM output file ", CASE);
  if (oid) SP_scan_expected_number_exit (str,         o, "Invalid structure of OOFEM output file");
  if (oid) SP_scan_expected_word_exit   (str,       "(", "Invalid structure of OOFEM output file", CASE);
  ;        SP_scan_expected_number_exit (str,         i, "Invalid structure of OOFEM output file");
  if (oid) SP_scan_expected_word_exit   (str,       ")", "Invalid structure of OOFEM output file", CASE);
}
void OOFEM_output_scan_GP (const char *&str, int i, bool old=false)
{
  char GP[7];
  if (old) sprintf (GP,   "%d", i);
  else     sprintf (GP, "1.%d", i);
  
  SP_scan_expected_word_exit   (str, "GP", "Invalid structure of OOFEM output file", CASE);
  SP_scan_expected_word_exit   (str,  GP , "Invalid structure of OOFEM output file", CASE);
  SP_scan_expected_word_exit   (str,  ":", "Invalid structure of OOFEM output file", CASE);
}
void OOFEM_output_scan_GP (FILE *stream, int i, bool old=false)
{
  char GP[7];
  if (old) sprintf (GP,   "%d", i);
  else     sprintf (GP, "1.%d", i);
  
  FP_scan_expected_word_exit   (stream, "GP", "Invalid structure of OOFEM output file", CASE);
  FP_scan_expected_word_exit   (stream,  GP , "Invalid structure of OOFEM output file", CASE);
  FP_scan_expected_word_exit   (stream,  ":", "Invalid structure of OOFEM output file", CASE);
}


//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS BEAM                          *** *** *** ***
//*  ********************************************************************************************
void Beam :: checkConsistency (void) const
{
  FElement :: checkConsistency();
  
  //* Detects beam with hinges/condensed rotation DOFs at both end point.
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_DxyzR___: break;
  case DPN_DxyzRxyz:
    if (this->elemAttribs()->has_fullhinge_ends())
      _errorr2("Rotation DOFs condensed at both ends of Beam - forbidden kinematic mechanism at elemnet %ld.", ID+1);
    break;
  default: errol;
  }
  
  // NA CO JE TOTO?
  // to bych mel spis kontrolovat na vyssi urovni, asi na Elemnetu, ale je to vubec nutne?
  // DOFsPerNode DPN = Pd->give_DOFspnod();
  //switch (this->elemAttribs()->give_dpn()) {
  //case DPN_DxyzR___:
  //case DPN_DxyzRxyz: if (DPN != DPN_DxyzR___  &&  DPN != DPN_DxyzRxyz) errol;  break;
  //default:  errol;
  //}
}

DOFsPerNode Beam :: give_DOFsPerNode_default (void) const
{
  switch (Pd->give_DOFspnod()) {
  case DPN_DxyzR___:  return DPN_DxyzR___;  break;
  case DPN_DxyzRxyz:  return DPN_DxyzRxyz;  break;
  default:            errol;
  }
  return DPN_Void;
}
SStype Beam :: give_SSType_default (void) const
{
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_DxyzR___:  return SST_truss;  break;
  case DPN_DxyzRxyz:  return SST_beam;   break;
  default:            errol;
  }
  return SST_Void;
}

// /// return type of element for OOFEM solver
// FiniteElementType Beam :: give_fe (Solver sol) const
// {
//   FiniteElementType feoDef, feoExpl;
// 
//   if      (elemAttribs()->give_rot() == false) feoDef = FET_Truss3d;
//   else                                      feoDef = FET_Beam3d;
// 
//   feoExpl = elemAttribs()->give_FET();
// 
//   if (feoExpl != FET_Void && feoExpl != feoDef) {
//     if      (feoExpl == FET_LiBeam3d    &&  feoDef == FET_Beam3d)  feoDef = FET_LiBeam3d;
//     else if (feoExpl == FET_LiBeam3dNL  &&  feoDef == FET_Beam3d)  feoDef = FET_LiBeam3dNL;
//     else errol;
//   }
// 
//   return feoDef;
// }

long Beam :: give_edge_nodes (const Point **&edgnodes) const
{
  if (elemAttribs()->give_ord() != -4) errol;
  edgnodes = new const Point*[2];

  edgnodes[0] = points[0];
  edgnodes[1] = points[1];
  return 2;
}

void Beam :: give_vector (VectoR *v) const
{
  if (points() != 2) errol;
  v->beP2P(points[0]->give_coords(), points[1]->give_coords());
}

void Beam :: read_output_OOFEM (FILE *stream, long step)
{
  const char *str;
  char LINE[1023];
  double *vald, l;
  Dmtrx *strain = new Dmtrx;
  Dmtrx *stress = new Dmtrx;
  
  IntPointSet ips = this->give_IPset_rslts(SOL_OOFEM);
  int NIP = IntPointSet_give_number_ips (ips);
  if (NIP != 1)  errol; // ip nacitat ve spravnem poradi, az nip != 1 !!!
  
  int ver = Pd->give_P_OOFEM_ver();
  
  FiniteElementTypeSet fets;
  FiniteElementType feto = FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_OOFEM, Pd->give_analgroup());
  
  switch (feto) {
  case FET_OOFEM_Truss3d:
    strain->resize_ignore_vals(NIP, 1);
    stress->resize_ignore_vals(NIP, 1);
    vald = new double[6];
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, ver >= 18);
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_GP (str, 1, ver < 18);
    SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
    if (! isZero(1e-10, vald[1]+vald[2]+vald[3]+vald[4]+vald[5]))  _errorr ("Invalid structure of OOFEM output file");
    (*strain)(0,0) = vald[0];
    
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
    if (! isZero(1e-10, vald[1]+vald[2]+vald[3]+vald[4]+vald[5]))  _errorr ("Invalid structure of OOFEM output file");
    (*stress)(0,0) = vald[0];

    delete [] vald;
    break;
  case FET_OOFEM_LiBeam3d:
    strain->resize_ignore_vals(NIP, 12);
    stress->resize_ignore_vals(NIP, 12);
    l = this->give_lav();

    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, ver >= 18);

    vald = strain->give_ptr2val(0, 0);
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_GP (str, 1, ver < 18);
    SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    if (! isZero(1e-10, vald[1]+vald[2]+vald[3]+vald[9]+vald[10]+vald[11]))  _errorr ("Invalid structure of OOFEM output file");
    vald[ 9] = l * vald[6];
    vald[10] = l * vald[7];
    vald[11] = l * vald[8];
    vald[ 6] = l * vald[0];
    vald[ 7] = l * vald[5];
    vald[ 8] = l * vald[4];
    vald[0] = vald[4] = vald[5] = 0.0;
    // ??
    vald[7] = vald[7] + 0.5 * l * vald[11];  // Vy += 0.5 * l * Mz
    vald[8] = vald[8] - 0.5 * l * vald[10];  // Vz -= 0.5 * l * My

    vald = stress->give_ptr2val(0, 0);
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    if (! isZero(1e-10, vald[1]+vald[2]+vald[3]+vald[9]+vald[10]+vald[11]))  _errorr ("Invalid structure of OOFEM output file");
    vald[1] = vald[5];
    vald[2] = vald[4];
    vald[3] = vald[ 9] = vald[6];
    vald[4] = vald[10] = vald[7];
    vald[5] = vald[11] = vald[8];
    vald[6] = vald[0];
    vald[7] = vald[1];
    vald[8] = vald[2];

    break;
  case FET_OOFEM_Beam3d:
    strain->resize_ignore_vals(NIP, 12);
    stress->resize_ignore_vals(NIP, 12);

    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit (str, "beam",    "Invalid structure of OOFEM output file", CASE);
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, false);
    
    // /// strain
    // vald = strain->give_ptr2val(0, 0);
    // fgets (LINE, 1023, stream);  str=LINE;
    // SP_scan_expected_word_exit (str, "local displacements", "Invalid structure of OOFEM output file", CASE);
    // if (!SP_scan_array (str, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    // /// stress
    // vald = stress->give_ptr2val(0, 0);
    // fgets (LINE, 1023, stream);  str=LINE;
    // SP_scan_expected_word_exit (str, "local end forces", "Invalid structure of OOFEM output file", CASE);
    // if (!SP_scan_array (str, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    // for (int i=0; i<6; i++)  vald[i] *= -1;
    
    vald = strain->give_ptr2val(0, 0);
    FP_scan_expected_word_exit (stream, "local displacements", "Invalid structure of OOFEM output file", CASE);
    if (!FP_scan_array (stream, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    
    vald = stress->give_ptr2val(0, 0);
    FP_scan_expected_word_exit (stream, "local end forces", "Invalid structure of OOFEM output file", CASE);
    if (!FP_scan_array (stream, 12, vald))  _errorr ("Invalid structure of OOFEM output file");
    for (int i=0; i<6; i++)  vald[i] *= -1;
    
    FP_skip_line(stream);
    
    break;
  default:  _errorr("unsupported type of beam element");
  }
  
  this->add_result(strain, step, RTE_global_strain);
  this->add_result(stress, step, RTE_global_stress);
}
void Beam :: read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt)
{
  const char *str;
  char LINE[1023], msg[60];
  
  sprintf (msg, "Invalid structure of SIFEL output file, element %ld", origid+1);
  
  IntPointSet ips = this->give_IPset_rslts(SOL_SIFEL);
  int NIP = IntPointSet_give_number_ips (ips);
  if (NIP != 1)  errol;  // toto je zde natvrdo, u tr v transportu je jeden IP,  // ip nacitat ve spravnem poradi, az nip != 1 !!!
  int NIPa = 1;          // ale tisknou se 4, protoze datel a nebo kaser a labus

  fgets (LINE, 1023, stream);  str=LINE;
  SP_scan_expected_word_exit   (str, "Element", msg, CASE);
  SP_scan_expected_number_exit (str,  origid+1, msg);
  SP_scan_expected_word_exit   (str, ", integration points", msg, CASE);
  long ip1, ip2;
  SP_scan_number (str, ip1);
  SP_scan_expected_word_exit   (str, "-", msg, CASE);
  SP_scan_number (str, ip2);
  if (ip2-ip1+1 != NIPa) _errorr(msg);

  int nstrcomp = 1;
  if (Pd->give_analgroup() == PAG_mechanics) {
    int j;
    Dmtrx *values = new Dmtrx(NIP, nstrcomp);

    // loop over NIP
    for (int i=0; i<NIP; i++) {
      fgets (LINE, 1023, stream);  str=LINE;

      char r[15];
      for (j=0; j<nstrcomp; j++) {
    if      (rt == RTE_global_strain)  sprintf(r, "eps_%d=", j+1);
    else if (rt == RTE_global_stress)  sprintf(r, "sig_%d=", j+1);
    else errol;
    SP_scan_expected_word_exit (str, r, msg, CASE);
    SP_scan_number (str, (*values)(i,j));
      }
    }

    FP_skip_line (stream, NIPa-NIP);

    this->add_result (values, step, rt);
  }
  else errol;
}

double Beam :: give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node)
{
  if (this->elemAttribs()->give_sst() != SST)  return this->fillupbyzero (data, SST);

  if ( Msh()->give_rslts_solver() == SOL_OOFEM ) {
    //if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
    if (node)  errol;

    if (SST == SST_truss) {
      data->resize_ignore_vals(1);

      if      (rvtype == RVTstrn) { check_rslts(step, RTE_global_strain);  data[0][0] = this->give_results_dm(step, RTE_global_strain)[0](0,0);  }
      else if (rvtype == RVTstrs) { check_rslts(step, RTE_global_stress);  data[0][0] = this->give_results_dm(step, RTE_global_stress)[0](0,0);  }
      else                        errol;
    }
    else if (SST == SST_beam) {
      data->resize_ignore_vals(9);

      const Dmtrx *s=NULL;

      if      (rvtype==RVTstrn) { check_rslts(step, RTE_global_strain);  s = this->give_results_dm(step, RTE_global_strain); }
      else if (rvtype==RVTstrs) { check_rslts(step, RTE_global_stress);  s = this->give_results_dm(step, RTE_global_stress); }
      else                      errol;

      if (rvtype==RVTstrn)  for (int i=0; i<6; i++)  data[0][i+3] =   s[0](0,i+6) - s[0](0,i)        ;
      else                  for (int i=0; i<6; i++)  data[0][i+3] = ( s[0](0,i+6) + s[0](0,i) ) / 2.0;

      data[0][0] = data[0][3];   // N
      data[0][1] = sqrt( data[0][4] * data[0][4]  +  data[0][5] * data[0][5] );   // V
      data[0][2] = sqrt( data[0][7] * data[0][7]  +  data[0][8] * data[0][8] );   // M
    }
    else errol;

    // bf==BF_Beam
    if (node)  _errorr("error");
    //bool first;
    //if (node == nodes[0])  first = true;
    //if (node == nodes[1])  first = false;
    //
    //switch (rvtype) {
    //case RVT_endDispl:
    //  if (!strain)  _errorr("no strain");
    //  if (first)  data->copy(this->give_results_dm(step, RTE_strain)->give_ptr2val()  , 6);
    //  else      { data->copy(this->give_results_dm(step, RTE_strain)->give_ptr2val()+6, 6);  data->tmsby(-1.0); }
    //  break;
    //case RVT_endForce:
    //  if (first)  data->copy(this->give_results_dm(step, RTE_stress)->give_ptr2val()  , 6);
    //  else      { data->copy(this->give_results_dm(step, RTE_stress)->give_ptr2val()+6, 6);  data->tmsby(-1.0); }
    //  break;
  }
  else if ( Msh()->give_rslts_solver() == SOL_SIFEL ) {
    //if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
    if (node)  errol;
    
    if (SST == SST_truss) {
      data->resize_ignore_vals(1);
      
      if      (rvtype == RVTstrn) { check_rslts(step, RTE_global_strain);  data[0][0] = this->give_results_dm(step, RTE_global_strain)[0](0,0);  }
      else if (rvtype == RVTstrs) { check_rslts(step, RTE_global_stress);  data[0][0] = this->give_results_dm(step, RTE_global_stress)[0](0,0);  }
      else                        errol;
    }
    else errol;
  }
  else errol;
  
  return this->give_GeomWeight1deg();
}

//double Beam :: give_normal_force (long step) const
//{
//  check_rslts(step);
//
//  Dmtrx *s = this->give_results_dm(step, RTE_stress);
//
//  if (bf == BF_Truss || bf == BF_LiBeam) {
//    return s[0](0,0);
//  }
//  else if (bf == BF_Beam) {
//    return ( s[0](0,6) - s[0](0,0) ) / 2.0;
//  }
//  else _errorr("unsupported type of beam element");
//
//  return 0.0;
//}
double SigmaEq (VectoR &s, VectoR &t)
{
  return 0.5 * sqrt ( 2.0 * ( (s.x-s.y)*(s.x-s.y) +
                  (s.y-s.z)*(s.y-s.z) +
                  (s.z-s.x)*(s.z-s.x) +
                  6.0 * (t.x*t.x + t.y*t.y + t.z*t.z)));
}

void Beam :: compute_maxSigmaEq ()
{
  long step = 0;
  double answer = 0.0;
  
  check_rslts(step, RTE_global_stress);
  
  if (this->elemAttribs()->give_sst() == SST_truss) {
    answer = this->give_results_dm(step, RTE_global_stress)[0](0,0);
    if (answer < 0)  answer *= -1;
  }
  else if (this->elemAttribs()->give_sst() == SST_beam) {
    if (this->elemAttribs()->give_cs()->give_loctype_or_type() == CST_Rectangle) {  // _or_type urizni
      double A, y, z, Iy, Iz;
      A =       this->elemAttribs()->give_cs()->give_area();
      y = 0.5 * this->elemAttribs()->give_cs()->give_width();
      z = 0.5 * this->elemAttribs()->give_cs()->give_height();
      Iy=       this->elemAttribs()->give_cs()->give_Iy();
      Iz=       this->elemAttribs()->give_cs()->give_Iz();
      
      VectoR s, t;
      double N, Vy, Vz, Mx, My, Mz, val;
      
      // 2  points == start & end
      for (int i=0; i<2; i++) {
    N  = fabs(this->give_results_dm(step, RTE_global_stress)[0](0,0+i*6) / A      );   //     ^ z       //
    Vy =      this->give_results_dm(step, RTE_global_stress)[0](0,1+i*6) / A * 1.5 ;   //  +--|--+      //
    Vz =      this->give_results_dm(step, RTE_global_stress)[0](0,2+i*6) / A * 1.5 ;   //  |  |  |       //
    Mx =      this->give_results_dm(step, RTE_global_stress)[0](0,3+i*6)           ;   //  3  4-----> y //
    My = fabs(this->give_results_dm(step, RTE_global_stress)[0](0,4+i*6) / Iy * z );   //  |     |       //
    Mz = fabs(this->give_results_dm(step, RTE_global_stress)[0](0,5+i*6) / Iz * y );   //  1--2--+       //
    
    // invariants
    s.y = s.z = 0.0;
    t.x = 0.0*Mx;   // Mx  docasne
    
    s.x = N + My + Mz;  t.y =  0;  t.z =  0;   val = SigmaEq(s, t);   if (answer < val)  answer = val;  // point 1
    s.x = N + My     ;  t.y =  0;  t.z = Vy;   val = SigmaEq(s, t);   if (answer < val)  answer = val;  // point 2
    s.x = N      + Mz;  t.y = Vz;  t.z =  0;   val = SigmaEq(s, t);   if (answer < val)  answer = val;  // point 3
    s.x = N          ;  t.y = Vz;  t.z = Vy;   val = SigmaEq(s, t);   if (answer < val)  answer = val;  // point 4
      }
    }
  }
  else errol;
  
  maxSigmaEq[0] = answer;
}

void Beam :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  FElement :: print_row (stream, fff, false);

  if (fff == FFF_OOFEM) {
    if (elemAttribs()->give_conDOFs()) {
      FiniteElementTypeSet fets;
      if (FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_OOFEM, Pd->give_analgroup()) != FET_OOFEM_Beam3d)  _errorr ("OOFEM supports hinges at Beam3d elements only");
      if (points() != 2)  errol;

      const GPA<PointDOFsCondense> *con = elemAttribs()->give_conDOFs();

      int i,j, n=0, ndofs=0;
      if      ((*con)[0]) ndofs = (*con)[0]->give_ndofs();
      else if ((*con)[1]) ndofs = (*con)[1]->give_ndofs();
      else errol;

      for (i=0; i<2; i++) { if (!(*con)[i]) continue; for (j=0; j<ndofs; j++) if ((*con)[i]->give_att(j)) n++; }
      fprintf (stream, "  dofsToCondense %d", n);
      for (i=0; i<2; i++) { if (!(*con)[i]) continue; for (j=0; j<ndofs; j++) if ((*con)[i]->give_att(j)) fprintf (stream, "  %d", i*ndofs+j+1); }
    }
  }
  else
    if (elemAttribs()->give_conDOFs()) errol;

  if (endline)  fprintf (stream,"\n");
}



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS TRIANGLE                        *** *** *** ***
//*  ********************************************************************************************
Triangle :: Triangle (const Quadrangle * src, bool firstr, bool firstdiag) :
  FElement (classComponentGeometry2Dtriangle, src->give_ID(), src->give_origID(), src->give_Geom(),1  ,3 ,3 ,1 , true, src->give_domain(), src->give_ID())
{
  if (Msh()->connectivity_is_assembled())  errol;
  
  this->mproperty = *(src->give_mproperty_ptr());
  if (src->give_delete_flag()) errol;
  
  const GPA<const Point> *nds = src->give_points();
  if (firstr) {
    points.assign(0, (*nds)[0]);
    points.assign(1, (*nds)[1]);
    if (firstdiag) points.assign(2, (*nds)[2]);
    else           points.assign(2, (*nds)[3]);
  }
  else {
    points.assign(0, (*nds)[2]);
    points.assign(1, (*nds)[3]);
    if (firstdiag) points.assign(2, (*nds)[0]);
    else           points.assign(2, (*nds)[1]);
  }
  
  if (src->give_domain()) errol;
  mdl_masterel = src->give_mdl_masterel();
  delete attributes;
  attributes = new FElemAttribs(this, (const FElemAttribs*)src->give_elemAttribs(), 'd');
  //attributes->set_IPset();
  if (elemAttribs()->give_prop())  errol;

  // this
  rotMg2l = NULL;
}

void Triangle :: checkConsistency (void) const
{
  FElement :: checkConsistency();
  
  FiniteElementTypeSet fets;
  this->elemAttribs()->give_FETS(&fets);
  
  switch (Pd->give_analgroup()) {
  case PAG_mechanics:
    ;    if (fets.dpn == DPN_Dxy_R___) { if (fets.sst != SST_plstrain  &&  fets.sst != SST_plstress)  errol; }
    else if (fets.dpn == DPN_Dxy_R__z) { if (fets.sst != SST_plstrain  &&  fets.sst != SST_plstress)  errol; }
    else if (fets.dpn == DPN_D__zRxy_) { if (fets.sst != SST_plate)                                   errol; }
    else if (fets.dpn == DPN_DxyzRxyz) { if (fets.sst != SST_3dshell)                                 errol; }
    else errol;
    break;
  case PAG_transport:
    if (fets.sst != SST_transp2d) errol;
    break;
  default: errol;
  }
}

DOFsPerNode Triangle :: give_DOFsPerNode_default (void) const
{
  switch (Pd->give_DOFspnod()) {
  case DPN_Dxy_R___:  return DPN_Dxy_R___;  break;
  case DPN_Dxy_R__z:  return DPN_Dxy_R__z;  break;
  case DPN_D__zRxy_:  return DPN_D__zRxy_;  break;
  case DPN_DxyzRxyz:  return DPN_DxyzRxyz;  break;
  default:            errol;
  }
  return DPN_Void;
}
SStype Triangle :: give_SSType_default (void) const
{
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_Dxy_R__z:  return SST_plstress;  break;
  case DPN_Dxy_R___:  return SST_plstress;  break;
  case DPN_D__zRxy_:  return SST_plate;     break;
  case DPN_DxyzRxyz:  return SST_3dshell;   break;
  default:            errol;
  }
  return SST_Void;
}

// FiniteElementType Triangle :: give_fe (Solver sol) const
// {
//   if (sol == SOL_OOFEM) {
//     switch (elemAttribs()->give_sst()) {
//     case SST_3d:           return FET_ShellTr;        break;
//     case SST_planestrain:  return FET_PlaneStrainTr;  break;
//     case SST_planestress:  return FET_PlaneStressTr;  break;
//     default:               errol;
//     }
//   }
//   if (sol == SOL_SIFEL)    return FET_PlaneTr;
//   
//   errol;  return FET_Void;
// } // "RerShell"; }


long Triangle :: give_edge_nodes (const Point **&edgnodes) const
{
  if (elemAttribs()->give_ord() != -4) _errorr2("ord is %d", elemAttribs()->give_ord());
  edgnodes = new const Point*[6];

  edgnodes[0] = points[0];   edgnodes[1] = points[1];
  edgnodes[2] = points[1];   edgnodes[3] = points[2];
  edgnodes[4] = points[2];   edgnodes[5] = points[0];
  return 2;
}
long Triangle :: give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const
{
  if (elemAttribs()->give_ord() != -4) errol; // musim vratit i pocet uzlu, to uz budu muset pres parametry
  facnodes = new const Point* [3];
  facedges = new const Edge* [3];

  facnodes[0] = points[0];  facnodes[ 1] = points[1];    facnodes[ 2] = points[2];
  facedges[0] = edges[0];   facedges[ 1] = edges[1];     facedges[ 2] = edges[2];
  return 3;
}

double Triangle :: give_quality (void) const
{
  if (elemAttribs()->give_ord() != -4)  return -1;

  // f = 4*root(3) = 6.9282032302755088
  // quality = f * area / (a*a + b*b + c*c);
  // return  6.9282032302755088 * this->give_lav() /
  //   ( this->give_edge(0)->give_lav() * this->give_edge(0)->give_lav() +
  //     this->give_edge(1)->give_lav() * this->give_edge(1)->give_lav() +
  //     this->give_edge(2)->give_lav() * this->give_edge(2)->give_lav()   ) ;

  double l1 = this->give_edge(0)->give_lav();
  double l2 = this->give_edge(1)->give_lav();
  double l3 = this->give_edge(2)->give_lav();

  return  6.9282032302755088 * this->give_lav() / ( l1*l1 + l2*l2 + l3*l3 );
}

void Triangle :: compute_rotMg2l (void)
{
  if (rotMg2l)  return;

  VectoR x, y, z;
  //
  x.beP2P( points[0]->give_coords(), points[1]->give_coords() );
  x.normalize();
  //
  y.beP2P( points[0]->give_coords(), points[2]->give_coords() );
  z.beVectProduct( &x, &y );
  z.normalize();
  //
  y.beVectProduct( &z, &x );
  y.normalize();

  //
  rotMg2l = new Dmtrx(3,3);
  //
  rotMg2l->copy_row(0,&x);
  rotMg2l->copy_row(1,&y);
  rotMg2l->copy_row(2,&z);
}

void Triangle :: rotate2local(Dvctr *data)
{
  //this->compute_rotMg2l();

  //data->tnsrRotWith(*rotMg2l);

  double a = 0.5 * atan ( 2*data[0][5] / (data[0][1] - data[0][0]) );

  if (fabs(a) > 1.e-6)
    data->tnsrRotAxisZangle (a);

  if (data[0][0] < data[0][1])
    data->tnsrRotAxisZangle (1.570796327);
}

void Triangle :: give_ip_coords_global (IntPointSet ips, int i, PoinT &coords) const
{
  PoinT l;  // local coords
  
  IPS_give_ip_coord_native (i, ips, l);
  
  coords.x = l[0] * points[0]->give_coord(0) + l[1] * points[1]->give_coord(0) + l[2] * points[2]->give_coord(0);
  coords.y = l[0] * points[0]->give_coord(1) + l[1] * points[1]->give_coord(1) + l[2] * points[2]->give_coord(1);
  coords.z = l[0] * points[0]->give_coord(2) + l[1] * points[1]->give_coord(2) + l[2] * points[2]->give_coord(2);
}

void Triangle :: read_input (const char *&str, femFileFormat fff)
{
  if (fff == FFF_SIFEL) {
    FiniteElementTypeSet fets;
    if (FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_SIFEL, Pd->give_analgroup()) == FET_SIFEL_planeelementlt) {
      int sst;
      SP_scan_number (str, sst);
      elemAttribs()->set_sst(SStype_i2e_SIFEL(sst));
    }
  }
  
  FElement :: read_input (str, fff);
}

void Triangle :: read_output_OOFEM (FILE *stream, long step)
{
//#ifdef DEBUG
  //if (! isZero(1e-10, vald[2]+vald[8]+vald[9]+vald[10]))  _errorr ("Invalid structure of OOFEM output file");
//#endif
  const char *str;
  char LINE[1023];
  int ipid;
  Dmtrx *strain = new Dmtrx;
  Dmtrx *stress = new Dmtrx;
  
  IntPointSet ips = this->give_IPset_rslts(SOL_OOFEM);
  int NIP = IntPointSet_give_number_ips (ips);
  int nl = this->elemAttribs()->give_cs()->give_nlayers();
  
  FiniteElementTypeSet fets;
  FiniteElementType feto = FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_OOFEM, Pd->give_analgroup());
  femFileFormat fff = FFF_OOFEM;
  
  switch (feto) {
  case FET_OOFEM_RerShell:
  case FET_OOFEM_Tr_Shell01: {
    bool rer = (feto == FET_OOFEM_RerShell) ? true : false;
    
    int ver = Pd->give_P_OOFEM_ver();
    
    //** specialni vypis ve lokalni verzi oofemu a v git verzi pro nip 1
    if (!rer && (ver == 20  ||  (ver >= 22  &&  NIP == 1))) {
      if (NIP != 1)  errol; // ip nacitat ve spravnem poradi, az nip != 1 !!!
      int i;
      
      // verze je specialni v tom, ze tam jsou krome globalnich i lokalni hodnoty
      Dmtrx *strain_loc = new Dmtrx;
      Dmtrx *stress_loc = new Dmtrx;
      
      strain->    resize_ignore_vals(NIP, 12);
      stress->    resize_ignore_vals(NIP, 12);
      strain_loc->resize_ignore_vals(NIP, 12);
      stress_loc->resize_ignore_vals(NIP, 12);
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, true);
      
      OOFEM_output_scan_GP (stream, 1);
      
      
      FP_scan_expected_word_exit (stream, "generalized strains", "Invalid structure of OOFEM output file", CASE);
      if (!FP_scan_array (stream, 12, strain->give_ptr2val()))  _errorr ("Invalid structure of OOFEM output file");
      for (i=6; i<12; i++)
    if (2 < i%6 && i%6 < 6)  (*strain)(0,i) /= 2.0;
      
      FP_scan_expected_word_exit (stream, "local", "Invalid structure of OOFEM output file", CASE);
      if (!FP_scan_array (stream, 12, strain_loc->give_ptr2val()))  _errorr ("Invalid structure of OOFEM output file");
      for (i=6; i<12; i++)
    if (2 < i%6 && i%6 < 6)  (*strain_loc)(0,i) /= 2.0;
      
# ifdef DEBUG
      // orotuje strain a porovna ho s strain_loc
      Dvctr vald(6);
      double l, auxd;
      this->compute_rotMg2l();
      
//       // D
//       vald.beCopyOf(strain->give_ptr2val(), 6);
//       vald.tnsrRotWith(*rotMg2l);
//       l = vald.give_lenght();
//       auxd = fabs(vald[2]);
//       if (l > 2e-8  &&  auxd > 1e-10  &&  !isZero(l*1e-8, auxd))
//  _errorr3 ("Invalid structure of OOFEM output file, l = %e, auxd = 0 = %e", l, auxd);  // 1e-4 is equal to precision of oofem strain/stress output
//       vald[2] = 0.0;
//       
//       for (int i=0; i<6; i++)
//  if (vald[i] > 0.0  &&  fabs(vald[i] - (*strain_loc)(0,i)) > 2e-10  &&  ! isZero(l*1.e-8, fabs((*strain_loc)(0,i))*1e-4, vald[i] - (*strain_loc)(0,i)))
//    _errorr4 ("Invalid structure of OOFEM output file, l = %e, %e  -  %e ", l, (*strain_loc)(0,i), auxd);  // 1e-4 is equal to precision of oofem strain/stress output
//       
//       // R
//       vald.beCopyOf(strain->give_ptr2val(0,6), 6);
//       vald.tnsrRotWith(*rotMg2l);
//       l = vald.give_lenght();
//       auxd = fabs(vald[2]) + fabs(vald[3]) + fabs(vald[4]);
//       if (l > 2e-8  &&  auxd > 1e-09  &&  !isZero(l*1e-8*3, auxd))
//  _errorr3 ("Invalid structure of OOFEM output file, l = %e, auxd = 0 = %e", l, auxd);  // 1e-4 is equal to precision of oofem strain/stress output
//       vald[2] = vald[3] = vald[4] = 0.0;
//       
//       for (int i=0; i<6; i++)
//  if (vald[i] > 0.0  &&  fabs(vald[i] - (*strain_loc)(0,i+6)) > 2e-10  &&  ! isZero(l*1.e-8, fabs((*strain_loc)(0,i+6))*1e-4, vald[i] - (*strain_loc)(0,i+6)))
//    _errorr4 ("Invalid structure of OOFEM output file, l = %e, %e  -  %e ", l, (*strain_loc)(0,i+6), auxd);  // 1e-4 is equal to precision of oofem strain/stress output
# endif
      
      //* STRESS
      
      //* global
      FP_scan_expected_word_exit (stream, "generalized stresses", "Invalid structure of OOFEM output file", CASE);
      if (!FP_scan_array (stream, 12, stress->give_ptr2val()))  _errorr ("Invalid structure of OOFEM output file");
      
      //* local
      FP_scan_expected_word_exit (stream, "local", "Invalid structure of OOFEM output file", CASE);
      if (!FP_scan_array (stream, 12, stress_loc->give_ptr2val()))  _errorr ("Invalid structure of OOFEM output file");
      
# ifdef DEBUG
      // orotuje strain a porovna ho s strain_loc
      
      // D
      vald.beCopyOf(stress->give_ptr2val(), 6);
      vald.tnsrRotWith(*rotMg2l);
      l = vald.give_lenght();
      auxd = fabs(vald[2]);
      if (l > 2e-8  &&  !isZero(l*1e-8, auxd)) //  &&  auxd > 1e-10
    _errorr3 ("Invalid structure of OOFEM output file, l = %e, auxd = 0 = %e", l, auxd);  // 1e-4 is equal to precision of oofem strain/stress output
      vald[2] = 0.0;
      
      for (int i=0; i<6; i++)
    if (vald[i] > 0.0  &&  fabs(vald[i] - (*stress_loc)(0,i)) > 2e-10  &&  ! isZero(l*1.e-8, fabs((*stress_loc)(0,i))*1e-4, vald[i] - (*stress_loc)(0,i)))
      _errorr4 ("Invalid structure of OOFEM output file, l = %e, %e  -  %e ", l, (*stress_loc)(0,i), auxd);  // 1e-4 is equal to precision of oofem strain/stress output
      
      // R
      vald.beCopyOf(stress->give_ptr2val(0,6), 6);
      vald.tnsrRotWith(*rotMg2l);
      l = vald.give_lenght();
      auxd = fabs(vald[2]) + fabs(vald[3]) + fabs(vald[4]);
      if (l > 2e-8  &&  auxd > 1e-09  &&  !isZero(l*1e-8*3, auxd))
    _errorr3 ("Invalid structure of OOFEM output file, l = %e, auxd = 0 = %e", l, auxd);  // 1e-4 is equal to precision of oofem strain/stress output
      vald[2] = vald[3] = vald[4] = 0.0;
      
      for (int i=0; i<6; i++)
    if (vald[i] > 0.0  &&  fabs(vald[i] - (*stress_loc)(0,i+6)) > 2e-10  &&  ! isZero(l*1.e-8, fabs((*stress_loc)(0,i+6))*1e-4, vald[i] - (*stress_loc)(0,i+6)))
      _errorr4 ("Invalid structure of OOFEM output file, l = %e, %e  -  %e ", l, (*stress_loc)(0,i+6), auxd);  // 1e-4 is equal to precision of oofem strain/stress output
# endif
      
      FP_skip_line(stream);
      
      this->add_result(strain_loc, step, RTE_localOO_strain);
      this->add_result(stress_loc, step, RTE_localOO_stress);
    }
    //** Strains and stresses in global coordinates.
    else {
      Dmtrx *lstrain;
      Dmtrx *lstress;
      
      if (nl) {
    lstrain = new Dmtrx;  lstrain->resize_ignore_vals(nl, 6);
    lstress = new Dmtrx;  lstress->resize_ignore_vals(nl, 6);
      }
      
      strain->resize_ignore_vals(NIP, 12);
      stress->resize_ignore_vals(NIP, 12);
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, !rer); // true);
      
      for (int i=0; i<NIP; i++) {
    ipid = ECN_convert_ips_triangle (false, i, fff, ips);
    
    //* STRAIN
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_GP (str, i+1, rer);
    SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 12, strain->give_ptr2val(ipid, 0)))  _errorr ("Invalid structure of OOFEM output file");
    
    for (int j=0; j<12; j++) // ????
      if (2 < j%6 && j%6 < 6)  (*strain)(ipid,j) /= 2.0;
    
    //* STRESS
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
    if (!SP_scan_array (str, 12, stress->give_ptr2val(ipid, 0)))  _errorr ("Invalid structure of OOFEM output file");
    
    //* LAYERS
    if (nl) {
      FP_scan_expected_line_exit(stream, "Layers report", "Invalid structure of OOFEM output file", CASE);
      FP_scan_expected_line_exit(stream, "{",             "Invalid structure of OOFEM output file", CASE);
      for (int j=0; j<nl; j++) {
        fgets (LINE, 1023, stream);  str=LINE;
        OOFEM_output_scan_GP (str, j+1);
        SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
        if (!SP_scan_array (str, 6, lstrain->give_ptr2val(j, 0)))  _errorr ("Invalid structure of OOFEM output file");
        
        fgets (LINE, 1023, stream);  str=LINE;
        SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
        if (!SP_scan_array (str, 6, lstress->give_ptr2val(j, 0)))  _errorr ("Invalid structure of OOFEM output file");
      }
      FP_scan_expected_line_exit(stream, "} end layers report",             "Invalid structure of OOFEM output file", CASE);
      
      this->add_result(lstrain, step, RTE_layered_strain);
      this->add_result(lstress, step, RTE_layered_stress);
    }
      }
    }
    
    break;
  }
  case FET_OOFEM_TrPlaneStress2d: {
    double vald[6];
    strain->resize_ignore_vals(NIP, 3);
    stress->resize_ignore_vals(NIP, 3);
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, true);
    
    for (int i=0; i<NIP; i++) {
      ipid = ECN_convert_ips_triangle (false, i, fff, ips);
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_GP (str, i+1);
      SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      (*strain)(ipid,0) = vald[0];
      (*strain)(ipid,1) = vald[1];
      (*strain)(ipid,2) = vald[5] / 2.0;
      
      fgets (LINE, 1023, stream);  str=LINE;
      SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      (*stress)(ipid,0) = vald[0];
      (*stress)(ipid,1) = vald[1];
      (*stress)(ipid,2) = vald[5];
    }
    
    if (elemAttribs()->give_mat()->isOOFEMplast()) {
      if (! FP_skip_expected_string (stream, "status { ", true))  _errorr("Invalid structure of OOFEM output file");
      if (FP_skip_expected_string (stream, " Yielding,  plastic strains", true)) {
    if (!this->give_results(step, RTE_plaststrain)) {
      this->add_result(new Dmtrx(NIP, 5), step, RTE_plaststrain);
      ((Dmtrx*)this->give_results_dm(step, RTE_plaststrain))->zero();
    }
    for (int j=0; j<4; j++)  fscanf (stream, "%lf", &(((Dmtrx*)this->give_results_dm(step, RTE_plaststrain))[0](0,j)) );
    //if (! FP_skip_expected_string (stream, ", strain space hardening vars", true))  _errorr3("Invalid structure of OOFEM output file, step %ld, element %ld", step+1, id+1);
    //fscanf (stream, "%lf", &((*plastthis->give_results_dm(step, RTE_strain))(i,4)) );
      }
      FP_skip_line (stream);
    }
    break;
  }
  case FET_OOFEM_TrPlaneStress2dRot: {
    // pro tento prvek se pouziva asi redukovana integrace, tj. je tam jedna sada navic
    // proto tady natvrdo kontroluju konkretni sadu nip
    if (NIP != 4) errol;
    
    double vald[6];
    strain->resize_ignore_vals(NIP, 3);
    stress->resize_ignore_vals(NIP, 3);
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, true);
    
    for (int i=0; i<NIP; i++) {
      ipid = ECN_convert_ips_triangle (false, i, fff, ips);
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_GP (str, i+1);
      SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      (*strain)(ipid,0) = vald[0];
      (*strain)(ipid,1) = vald[1];
      (*strain)(ipid,2) = vald[5] / 2.0;
      
      fgets (LINE, 1023, stream);  str=LINE;
      SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      (*stress)(ipid,0) = vald[0];
      (*stress)(ipid,1) = vald[1];
      (*stress)(ipid,2) = vald[5];
    }
    
    // preskakuje ip redukovany
    FP_skip_line (stream, 2);
    
    if (elemAttribs()->give_mat()->isOOFEMplast())  errol;
    break;
  }
  default:  _errorr("unsupported type of triangle element");
  }
  
  this->add_result(strain, step, RTE_global_strain);
  this->add_result(stress, step, RTE_global_stress);
}
void Triangle :: read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt)
{
  const char *str;
  char LINE[1023], msg[60];

  sprintf (msg, "Invalid structure of SIFEL output file, element %ld", origid+1);
  
  IntPointSet ips = this->give_IPset_rslts(SOL_SIFEL);
  int NIP = IntPointSet_give_number_ips (ips);
  if (NIP != 1)  errol;  // toto je zde natvrdo, u tr v transportu je jeden IP, // ip nacitat ve spravnem poradi, az nip != 1 !!!
  int NIPa = 4;          // ale tisknou se 4, protoze datel a nebo kaser a labus

  fgets (LINE, 1023, stream);  str=LINE;
  SP_scan_expected_word_exit   (str, "Element", msg, CASE);
  SP_scan_expected_number_exit (str,  origid+1, msg);
  SP_scan_expected_word_exit   (str, ", integration points", msg, CASE);
  long ip1, ip2;
  SP_scan_number (str, ip1);
  SP_scan_expected_word_exit   (str, "-", msg, CASE);
  SP_scan_number (str, ip2);
  if (ip2-ip1+1 != NIPa) _errorr(msg);

  if (Pd->give_analgroup() == PAG_transport) {
    int dim, nmed, j, k;
    dim = this->give_dimension();
    nmed = Pd->give_global_nDOFs();

    Dmtrx *values = new Dmtrx(NIP, dim*nmed);

    // loop over NIP
    fgets (LINE, 1023, stream);  str=LINE;

    char r[15];
    for (j=0; j<nmed; j++)
      for (k=0; k<dim; k++) {
    if      (rt == RTE_gradient)  sprintf(r, "grad_%d_%d=", j+1, k+1);
    else if (rt == RTE_fluxes)    sprintf(r, "flux_%d_%d=", j+1, k+1);
    else errol;
    SP_scan_expected_word_exit (str, r, msg, CASE);
    SP_scan_number (str, (*values)(0,j*dim+k));
      }

    FP_skip_line (stream, NIPa-NIP);

    this->add_result (values, step, rt);
  }
  else errol;
}

double Triangle :: give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node)
{
  if (this->elemAttribs()->give_sst() != SST)  return this->fillupbyzero (data, SST);
  
  if ( Msh()->give_rslts_solver() == SOL_OOFEM ) {
    if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
    
    if (SST == SST_3dshell) {
      switch (rvtype) {
      case RVTstrn_displ:             check_rslts(step, RTE_global_strain);   data->beCopyOf(this->give_results_dm(step, RTE_global_strain) ->give_ptr2val()  , 6);  break;
      case RVTstrn_rotat:             check_rslts(step, RTE_global_strain);   data->beCopyOf(this->give_results_dm(step, RTE_global_strain) ->give_ptr2val()+6, 6);  break;
      case RVTstrs_displ:             check_rslts(step, RTE_global_stress);   data->beCopyOf(this->give_results_dm(step, RTE_global_stress) ->give_ptr2val()  , 6);  break;
      case RVTstrs_rotat:             check_rslts(step, RTE_global_stress);   data->beCopyOf(this->give_results_dm(step, RTE_global_stress) ->give_ptr2val()+6, 6);  break;
      case RVTstrn_displ_local_main:  check_rslts(step, RTE_localOO_strain);  data->beCopyOf(this->give_results_dm(step, RTE_localOO_strain)->give_ptr2val()  , 6);  this->rotate2local(data);  break;
      case RVTstrn_rotat_local_main:  check_rslts(step, RTE_localOO_strain);  data->beCopyOf(this->give_results_dm(step, RTE_localOO_strain)->give_ptr2val()+6, 6);  this->rotate2local(data);  break;
      case RVTstrs_displ_local_main:  check_rslts(step, RTE_localOO_stress);  data->beCopyOf(this->give_results_dm(step, RTE_localOO_stress)->give_ptr2val()  , 6);  this->rotate2local(data);  break;
      case RVTstrs_rotat_local_main:  check_rslts(step, RTE_localOO_stress);  data->beCopyOf(this->give_results_dm(step, RTE_localOO_stress)->give_ptr2val()+6, 6);  this->rotate2local(data);  break;
      default:             errol;
      }
      //
    }
    else if (SST == SST_plstrain  ||  SST == SST_plstress) {
      if      (rvtype == RVTstrn_displ) { check_rslts(step, RTE_global_strain);  data->beCopyOf(this->give_results_dm(step, RTE_global_strain)->give_ptr2val(), 3); }
      else if (rvtype == RVTstrs_displ) { check_rslts(step, RTE_global_stress);  data->beCopyOf(this->give_results_dm(step, RTE_global_stress)->give_ptr2val(), 3); }
      else errol;
    }
    else errol;
  }
  else if ( Msh()->give_rslts_solver() == SOL_SIFEL ) {
    errol;
    //   if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
  //   if (rvtype != RVTstrs_grad) errol;
  // 
  //   if (SST == SST_transp2d) {
  //     int nmed = Pd->give_global_nDOFs();    if (nmed != 1  &&  nmed != 2)  errol;
  //     data->resize_ignore_vals(3*nmed);  data->zero();
  //     if (nmed>0) { (*data)[0] = this->give_results_dm(step, RTE_gradient)->give_ptr2val()[0];
  //                (*data)[1] = this->give_results_dm(step, RTE_gradient)->give_ptr2val()[1]; }
  //     if (nmed>1) { (*data)[3] = this->give_results_dm(step, RTE_gradient)->give_ptr2val()[2];
  //                (*data)[4] = this->give_results_dm(step, RTE_gradient)->give_ptr2val()[3]; }
  //   }
  //   else errol;
  //   
  }
  else errol;
  
  return this->give_GeomWeight1deg();
}

void Triangle :: compute_maxSigmaEq ()
{
  long step = 0;
  double answer = 0.0;
  
  if (this->elemAttribs()->give_sst() == SST_3dshell) {
    switch (this->elemAttribs()->give_cs()->give_loctype_or_type()) {
    case CST_2d: {
      check_rslts(step, RTE_localOO_stress);
      
      double h;
      h = this->elemAttribs()->give_cs()->give_thickness();
      
      VectoR s, t;
      double nx, ny, vyz, vxz, vxy, mx, my, mxy, val;
      
      nx  = this->give_results_dm(step, RTE_localOO_stress)[0](0, 0)       / h;
      ny  = this->give_results_dm(step, RTE_localOO_stress)[0](0, 1)       / h;
      
      vyz = this->give_results_dm(step, RTE_localOO_stress)[0](0, 3) * 1.5 / h;
      vxz = this->give_results_dm(step, RTE_localOO_stress)[0](0, 4) * 1.5 / h;
      
      vxy = this->give_results_dm(step, RTE_localOO_stress)[0](0, 5)       / h;
      
      mx  = this->give_results_dm(step, RTE_localOO_stress)[0](0, 6) * 6.0 / h / h;
      my  = this->give_results_dm(step, RTE_localOO_stress)[0](0, 7) * 6.0 / h / h;
      mxy = this->give_results_dm(step, RTE_localOO_stress)[0](0,11) * 6.0 / h / h;
  
      s.x = nx + mx;    t.x = 0;          
      s.y = ny + my;    t.y = 0;          
      s.z = 0;      t.z = vxy + mxy;
      
      val = SigmaEq(s, t);
      if (answer < val)  answer = val;
  
      s.x = nx - mx;    t.x = 0;          
      s.y = ny - my;    t.y = 0;          
      s.z = 0;          t.z = vxy - mxy;
      
      val = SigmaEq(s, t);
      if (answer < val)  answer = val;
  
      s.x = nx;         t.x = vyz;
      s.y = ny;         t.y = vxz;
      s.z = 0;          t.z = vxy;
      
      val = SigmaEq(s, t);
      if (answer < val)  answer = val;
      
      maxSigmaEq[0] = answer;
      break;
    }
    case CST_LayeredCS: {
      check_rslts(step, RTE_layered_stress);
      VectoR s, t;
      Dmtrx strs = this->give_results_dm(step, RTE_layered_stress);
      int nl = this->elemAttribs()->give_cs()->give_nlayers();
      
      for (int i=0; i<nl; i++) {
    s.x = strs(i,0);
    s.y = strs(i,1);
    s.z = strs(i,2);
    t.x = strs(i,3);
    t.y = strs(i,4);
    t.z = strs(i,5);
    
    maxSigmaEq[i] = SigmaEq(s, t);
      }
      break;
    }
    default: errol;
    }
  }
  else if (this->elemAttribs()->give_sst() == SST_plstress) {
    maxSigmaEq[0] = 0.0;
  }
  else errol;
}



//*  ********************************************************************************************
//*  *** *** *** ***                      CLASS QUADRANGLE                      *** *** *** ***
//*  ********************************************************************************************
void Quadrangle :: checkConsistency (void) const
{
  FElement :: checkConsistency();
  
  FiniteElementTypeSet fets;
  this->elemAttribs()->give_FETS(&fets);
  
  switch (Pd->give_analgroup()) {
  case PAG_mechanics:
    ;    if (fets.dpn == DPN_Dxy_R___) { if (fets.sst != SST_plstrain  &&  fets.sst != SST_plstress)  errol; }
    else if (fets.dpn == DPN_D__zRxy_) { if (fets.sst != SST_plate)                                    errol; }
    else if (fets.dpn == DPN_DxyzRxyz) { if (fets.sst != SST_3dshell)                                  errol; }
    else errol;  break;
  case PAG_transport:
    if (fets.sst != SST_transp2d) errol;  break;
  default: errol;
  }
}

DOFsPerNode Quadrangle :: give_DOFsPerNode_default (void) const
{
  switch (Pd->give_DOFspnod()) {
  case DPN_Dxy_R___:  return DPN_Dxy_R___;  break;
  case DPN_D__zRxy_:  return DPN_D__zRxy_;  break;
  case DPN_DxyzRxyz:  return DPN_DxyzRxyz;  break;
  default:            errol;
  }
  return DPN_Void;
}
SStype Quadrangle :: give_SSType_default (void) const
{
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_Dxy_R___:  return SST_plstress;  break;
  case DPN_D__zRxy_:  return SST_plate;     break;
  case DPN_DxyzRxyz:  return SST_3dshell;   break;
  default:            errol;
  }
  return SST_Void;
}




// FiniteElementType Quadrangle :: give_fe (Solver sol) const
// {
//   if (sol == SOL_OOFEM) {
//     if (elemAttribs()->give_sst() == SST_planestrain)  return FET_PlaneStrainQuad;
//     if (elemAttribs()->give_sst() == SST_planestress)  return FET_PlaneStressQuad;
//   }
//   if (sol == SOL_SIFEL)                             return FET_PlaneQuad;
// 
//   errol; return FET_Void;
// }


long Quadrangle :: give_edge_nodes (const Point **&edgnodes) const
{
  if (elemAttribs()->give_ord() != -4) errol;
  edgnodes = new const Point*[8];

  edgnodes[0] = points[0];   edgnodes[1] = points[1];
  edgnodes[2] = points[1];   edgnodes[3] = points[2];
  edgnodes[4] = points[2];   edgnodes[5] = points[3];
  edgnodes[6] = points[3];   edgnodes[7] = points[0];
  return 2;
}
long Quadrangle :: give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const
{
  if (elemAttribs()->give_ord() != -4) errol; // musim vratit i pocet uzlu, to uz budu muset pres parametry
  facnodes = new const Point* [4];
  facedges = new const Edge* [4];

  facnodes[0] = points[0];  facnodes[1] = points[1];  facnodes[2] = points[2];  facnodes[3] = points[3];
  facedges[0] = edges[0];   facedges[1] = edges[1];   facedges[2] = edges[2];   facedges[3] = edges[3];
  return 4;
}

bool Quadrangle :: is_first_diag_short (void)
{
  double d1 = points[0]->give_coords()->dist2_to(points[2]->give_coords());
  double d2 = points[1]->give_coords()->dist2_to(points[3]->give_coords());
  if (d1 < d2) return true;
  return false;
}

void Quadrangle :: read_input (const char *&str, femFileFormat fff)
{
  if (fff == FFF_SIFEL) {
    FiniteElementTypeSet fets;
    if (FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_SIFEL, Pd->give_analgroup()) == FET_SIFEL_planeelementlq) {
      int sst;
      SP_scan_number (str, sst);
      elemAttribs()->set_sst(SStype_i2e_SIFEL(sst));
    }
  }
  
  FElement :: read_input (str, fff);
}

void Quadrangle :: read_output_OOFEM (FILE *stream, long step)
{
  const char *str;
  char LINE[1023];
  Dmtrx *strain = new Dmtrx;
  Dmtrx *stress = new Dmtrx;
  
  IntPointSet ips = this->give_IPset_rslts(SOL_OOFEM);
  int NIP = IntPointSet_give_number_ips (ips);
  if (NIP != 4)  errol;
  
  FiniteElementTypeSet fets;
  this->elemAttribs()->give_FETS(&fets);
  FiniteElementType feto = FETSet_set2e (&fets, SOL_OOFEM, Pd->give_analgroup());
  
  switch (feto) {
  case FET_OOFEM_PlaneStress2d:
  case FET_OOFEM_Quad1PlaneStrain: {
    int i, j, ipid;
    double vald[6];
    strain->resize_ignore_vals(NIP, 3);
    stress->resize_ignore_vals(NIP, 3);
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, false);
    
    for (i=0; i<NIP; i++) {
      ipid = ECN_convert_ips_quadrangle (false, i, FFF_OOFEM, ips);   // id of ip; mapping oofem ipid to midas ipid
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_GP (str, i+1);
      SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      (*strain)(ipid,0) = vald[0];
      (*strain)(ipid,1) = vald[1];
      (*strain)(ipid,2) = vald[5] / 2.0;
      
      fgets (LINE, 1023, stream);  str=LINE;
      SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, vald))  _errorr ("Invalid structure of OOFEM output file");
#ifdef DEBUG
      if (fets.sst == SST_plstrain) vald[2] = 0.0;
      if (! isZero(1e-10, vald[2]+vald[3]+vald[4]))  _errorr ("Invalid structure of OOFEM output file");
#endif
      if (1) {
    (*stress)(ipid,0) = vald[0];
    (*stress)(ipid,1) = vald[1];
    (*stress)(ipid,2) = vald[5];
      }
      else {
    (*stress)(ipid,0) = vald[0] * elemAttribs()->give_cs()->give_thickness();
    (*stress)(ipid,1) = vald[1] * elemAttribs()->give_cs()->give_thickness();
    (*stress)(ipid,2) = vald[5] * elemAttribs()->give_cs()->give_thickness();
      }
      
      if (elemAttribs()->give_mat()->isOOFEMplast()) {
    if (! FP_skip_expected_string (stream, "status { ", true))  _errorr("Invalid structure of OOFEM output file");
    if (FP_skip_expected_string (stream, " Yielding,  plastic strains", true)) {
      if (!this->give_results(step, RTE_plaststrain)) {
        this->add_result(new Dmtrx(NIP, 5), step, RTE_plaststrain);
        ((Dmtrx*)this->give_results_dm(step, RTE_plaststrain))->zero();
      }
      for (j=0; j<4; j++)  fscanf (stream, "%lf", &(((Dmtrx*)this->give_results_dm(step, RTE_plaststrain))[0](ipid,j)) );
      // if (! FP_skip_expected_string (stream, ", strain space hardening vars", true))  _errorr3("Invalid structure of OOFEM output file, step %ld, element %ld", step+1, ID+1);
      // ;                    fscanf (stream, "%lf", &(((Dmtrx*)this->give_results_dm(step, RTE_plaststrain))[0](ipid,4)) );
    }
    FP_skip_line (stream);
      }
    }
    break;
  }
  default:  _errorr("unsupported type of brick element");
  }
  
  this->add_result(strain, step, RTE_global_strain);
  this->add_result(stress, step, RTE_global_stress);
}

double Quadrangle :: give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node)
{
  if (this->elemAttribs()->give_sst() != SST)  return this->fillupbyzero (data, SST);
  
  if ( Msh()->give_rslts_solver() == SOL_OOFEM ) {
    if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
    
    if (SST == SST_plstrain  ||  SST == SST_plstress) {
      if      (rvtype == RVTstrn_displ) { check_rslts(step, RTE_global_strain);  data->be_mean_of(this->give_results_dm(step, RTE_global_strain)); }
      else if (rvtype == RVTstrs_displ) { check_rslts(step, RTE_global_stress);  data->be_mean_of(this->give_results_dm(step, RTE_global_stress)); }
      else errol;
    }
    else errol;
  }
  else errol;
  
  if (type != 'g')  errol;
  //this->rotate2local(data);
  
  return this->give_GeomWeight1deg();
}

void Quadrangle :: give_ip_coords_global (IntPointSet ips, int i, PoinT &coords) const
{
  PoinT ncoords;  // local coords
  
  IPS_give_ip_coord_native (i, ips, ncoords);
  
  double l[4];
  
  l[0] = 0.25 * (1 + ncoords[0]) * (1 + ncoords[1]);
  l[1] = 0.25 * (1 - ncoords[0]) * (1 + ncoords[1]);
  l[2] = 0.25 * (1 - ncoords[0]) * (1 - ncoords[1]);
  l[3] = 0.25 * (1 + ncoords[0]) * (1 - ncoords[1]);
  
  coords.x = l[0] * points[0]->give_coord(0) + l[1] * points[1]->give_coord(0) + l[2] * points[2]->give_coord(0) + l[3] * points[3]->give_coord(0);
  coords.y = l[0] * points[0]->give_coord(1) + l[1] * points[1]->give_coord(1) + l[2] * points[2]->give_coord(1) + l[3] * points[3]->give_coord(1);
  coords.z = 0.0;
}


//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS TETRA                          *** *** *** ***
//*  ********************************************************************************************
void Tetra :: checkConsistency (void) const
{
  FElement :: checkConsistency();
  
  DOFsPerNode dpn = this->elemAttribs()->give_dpn();
  
  if (dpn != DPN_DxyzR___)  errol;
}

DOFsPerNode Tetra :: give_DOFsPerNode_default (void) const
{
  switch (Pd->give_DOFspnod()) {
  case DPN_DxyzR___:  return DPN_DxyzR___;  break;
  default:            errol;
  }
  return DPN_Void;
}
SStype Tetra :: give_SSType_default (void) const
{
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_DxyzR___:  return SST_3d;  break;
  default:            errol;
  }
  return SST_Void;
}


// ///
// FiniteElementType Tetra :: give_fe (Solver sol) const
// {
//   if (elemAttribs()->give_ord() != -4 || elemAttribs()->give_rot())  errol;
//   if (elemAttribs()->give_sst() == SST_3d)  return FET_TetraHedra;
//   else errol;  return FET_Void;
// }

void Tetra :: print_row_VTX (char *str) const
{
  sprintf (str, " %ld %ld %ld %ld",
       ((Node*)points[0])->give_lid_id(0),
       ((Node*)points[1])->give_lid_id(0),
       ((Node*)points[2])->give_lid_id(0),
       ((Node*)points[3])->give_lid_id(0));
}

long Tetra :: give_edge_nodes (const Point **&edgnodes) const
{
  if (elemAttribs()->give_ord() != -4) errol;
  edgnodes = new const Point*[12];

  edgnodes[ 0] = points[0];       edgnodes[ 1] = points[1];
  edgnodes[ 2] = points[1];       edgnodes[ 3] = points[2];
  edgnodes[ 4] = points[2];       edgnodes[ 5] = points[0];

  edgnodes[ 6] = points[0];       edgnodes[ 7] = points[3];
  edgnodes[ 8] = points[1];       edgnodes[ 9] = points[3];
  edgnodes[10] = points[2];       edgnodes[11] = points[3];

  return 2;
}
long Tetra :: give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const
{
  if (elemAttribs()->give_ord() != -4) errol; // musim vratit i pocet uzlu, to uz budu muset pres parametry
  facnodes = new const Point* [12];
  facedges = new const Edge*  [12];

  facnodes[ 0] = points[0];   facnodes[ 1] = points[1];   facnodes[ 2] = points[2];
  facnodes[ 3] = points[0];   facnodes[ 4] = points[1];   facnodes[ 5] = points[3];
  facnodes[ 6] = points[1];   facnodes[ 7] = points[2];   facnodes[ 8] = points[3];
  facnodes[ 9] = points[2];   facnodes[10] = points[0];   facnodes[11] = points[3];

  facedges[ 0] = edges[0];    facedges[ 1] = edges[1];    facedges[ 2] = edges[2];
  facedges[ 3] = edges[0];    facedges[ 4] = edges[4];    facedges[ 5] = edges[3];
  facedges[ 6] = edges[1];    facedges[ 7] = edges[5];    facedges[ 8] = edges[4];
  facedges[ 9] = edges[2];    facedges[10] = edges[3];    facedges[11] = edges[5];

  return 3;
}


//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS BRICK                          *** *** *** ***
//*  ********************************************************************************************
void Brick :: checkConsistency (void) const
{
  FElement :: checkConsistency();
  
  if (Pd->give_fulldata()) {
    DOFsPerNode dpn = this->elemAttribs()->give_dpn();
    
    if (dpn != DPN_DxyzR___)  errol;
  }
}

DOFsPerNode Brick :: give_DOFsPerNode_default (void) const
{
  switch (Pd->give_DOFspnod()) {
  case DPN_DxyzR___:  return DPN_DxyzR___;  break;
  default:            errol;
  }
  return DPN_Void;
}
SStype Brick :: give_SSType_default (void) const
{
  switch (this->elemAttribs()->give_dpn()) {
  case DPN_DxyzR___:  return SST_3d;  break;
  default:            errol;
  }
  return SST_Void;
}


// ///
// FiniteElementType Brick :: give_fe (Solver sol) const
// {
//   if (elemAttribs()->give_ord() != -4 || elemAttribs()->give_rot())  errol;
//   if (elemAttribs()->give_sst() == SST_3d)  return FET_HexaHedra;
//   else errol;  return FET_Void;  //return FET_QSpace; //return FET_LSpaceRot;
// }

long Brick :: give_edge_nodes (const Point **&edgnodes) const
{
  if (elemAttribs()->give_ord() != -4) errol;
  edgnodes = new const Point*[24];

  edgnodes[ 0] = points[0];       edgnodes[ 1] = points[1];
  edgnodes[ 2] = points[1];       edgnodes[ 3] = points[2];
  edgnodes[ 4] = points[2];       edgnodes[ 5] = points[3];
  edgnodes[ 6] = points[3];       edgnodes[ 7] = points[0];

  edgnodes[ 8] = points[4];       edgnodes[ 9] = points[5];
  edgnodes[10] = points[5];       edgnodes[11] = points[6];
  edgnodes[12] = points[6];       edgnodes[13] = points[7];
  edgnodes[14] = points[7];       edgnodes[15] = points[4];

  edgnodes[16] = points[0];       edgnodes[17] = points[4];
  edgnodes[18] = points[1];       edgnodes[19] = points[5];
  edgnodes[20] = points[2];       edgnodes[21] = points[6];
  edgnodes[22] = points[3];       edgnodes[23] = points[7];
  return 2;
}
long Brick :: give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const
{
  if (elemAttribs()->give_ord() != -4) errol; // musim vratit i pocet uzlu, to uz budu muset pres parametry
  facnodes = new const Point* [24];
  facedges = new const Edge* [24];

  facnodes[ 0] = points[0];   facnodes[ 1] = points[1];     facnodes[ 2] = points[2];   facnodes[ 3] = points[3];
  facnodes[ 4] = points[4];   facnodes[ 5] = points[5];     facnodes[ 6] = points[6];   facnodes[ 7] = points[7];
  facnodes[ 8] = points[0];   facnodes[ 9] = points[1];     facnodes[10] = points[5];   facnodes[11] = points[4];
  facnodes[12] = points[1];   facnodes[13] = points[2];     facnodes[14] = points[6];   facnodes[15] = points[5];
  facnodes[16] = points[2];   facnodes[17] = points[3];     facnodes[18] = points[7];   facnodes[19] = points[6];
  facnodes[20] = points[3];   facnodes[21] = points[0];     facnodes[22] = points[4];   facnodes[23] = points[7];

  facedges[ 0] = edges[0];   facedges[ 1] = edges[ 1];    facedges[ 2] = edges[2];   facedges[ 3] = edges[ 3];
  facedges[ 4] = edges[4];   facedges[ 5] = edges[ 5];    facedges[ 6] = edges[6];   facedges[ 7] = edges[ 7];
  facedges[ 8] = edges[0];   facedges[ 9] = edges[ 9];    facedges[10] = edges[4];   facedges[11] = edges[ 8];
  facedges[12] = edges[1];   facedges[13] = edges[10];    facedges[14] = edges[5];   facedges[15] = edges[ 9];
  facedges[16] = edges[2];   facedges[17] = edges[11];    facedges[18] = edges[6];   facedges[19] = edges[10];
  facedges[20] = edges[3];   facedges[21] = edges[ 8];    facedges[22] = edges[7];   facedges[23] = edges[11];
  return 4;
}

void Brick :: read_output_OOFEM (FILE *stream, long step)
{
  const char *str;
  char LINE[1023];
  Dmtrx *strain = new Dmtrx;
  Dmtrx *stress = new Dmtrx;
  
  IntPointSet ips = this->give_IPset_rslts(SOL_OOFEM);
  int NIP = IntPointSet_give_number_ips (ips);
  if (NIP != 8)  errol;
  
  FiniteElementTypeSet fets;
  FiniteElementType feto = FETSet_set2e (this->elemAttribs()->give_FETS(&fets), SOL_OOFEM, Pd->give_analgroup());
  
  switch (feto) {
  case FET_OOFEM_LSpace: {
    int i, ipid;
    strain->resize_ignore_vals(NIP, 6);
    stress->resize_ignore_vals(NIP, 6);
    
    int ver = Pd->give_P_OOFEM_ver();
    
    fgets (LINE, 1023, stream);  str=LINE;
    OOFEM_output_scan_elem_head (str, origid+1, this->ID+1, ver >= 18);
    
    for (i=0; i<NIP; i++) {
      ipid = ECN_convert_ips_brick (false, i, FFF_OOFEM, ips);   // id of ip; mapping oofem ipid to midas ipid
      
      fgets (LINE, 1023, stream);  str=LINE;
      OOFEM_output_scan_GP (str, i+1, ver < 18);
      SP_scan_expected_word_exit (str, "strains", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, strain->give_ptr2val(ipid, 0)))  _errorr ("Invalid structure of OOFEM output file");
      strain->give_ptr2val(ipid, 0)[3] = strain->give_ptr2val(ipid, 0)[3] / 2.0;
      strain->give_ptr2val(ipid, 0)[4] = strain->give_ptr2val(ipid, 0)[4] / 2.0;
      strain->give_ptr2val(ipid, 0)[5] = strain->give_ptr2val(ipid, 0)[5] / 2.0;
      
      fgets (LINE, 1023, stream);  str=LINE;
      SP_scan_expected_word_exit (str, "stresses", "Invalid structure of OOFEM output file", CASE);
      if (!SP_scan_array (str, 6, stress->give_ptr2val(ipid, 0)))  _errorr ("Invalid structure of OOFEM output file");
    }
    break;
  }
  default:  _errorr("unsupported type of brick element");
  }
  
  this->add_result(strain, step, RTE_global_strain);
  this->add_result(stress, step, RTE_global_stress);
}

void Brick :: read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt)
{
  const char *str;
  char LINE[1023], msg[60];
  
  sprintf (msg, "Invalid structure of SIFEL output file, element %ld", origid+1);
  
  int NIP = this->give_rslt_NIP(SOL_SIFEL);
  if (NIP != 8)  errol;  // toto je zde natvrdo
  int NIPa = 8;          // ...
  
  fgets (LINE, 1023, stream);  str=LINE;
  SP_scan_expected_word_exit   (str, "Element", msg, CASE);
  SP_scan_expected_number_exit (str,  origid+1, msg);
  SP_scan_expected_word_exit   (str, ", integration points", msg, CASE);
  long ip1, ip2;
  SP_scan_number (str, ip1);
  SP_scan_expected_word_exit   (str, "-", msg, CASE);
  SP_scan_number (str, ip2);
  if (ip2-ip1+1 != NIPa) _errorr(msg);
  
  int nstrcomp = 6;
  if (Pd->give_analgroup() == PAG_mechanics) {
    int j;
    Dmtrx *values = new Dmtrx(NIP, nstrcomp);
    
    // loop over NIP
    for (int i=0; i<NIP; i++) {
      // ipid = ECN_convert_ips_quadrangle (false, i, FFF_OOFEM);   // id of ip; mapping oofem ipid to midas ipid
      
      fgets (LINE, 1023, stream);  str=LINE;
      
      char r[15];
      for (j=0; j<nstrcomp; j++) {
    if      (rt == RTE_global_strain)  sprintf(r, "eps_%d=", j+1);
    else if (rt == RTE_global_stress)  sprintf(r, "sig_%d=", j+1);
    else errol;
    SP_scan_expected_word_exit (str, r, msg, CASE);
    SP_scan_number (str, (*values)(i,j));
      }
    }
    
    FP_skip_line (stream, NIPa-NIP);
    
    this->add_result (values, step, rt);
  }
  else errol;
}

void Brick :: read_input (const char *&str, femFileFormat fff)
{
  FElement :: read_input (str, fff);
}

double Brick :: give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node)
{
  if (this->elemAttribs()->give_sst() != SST)  return this->fillupbyzero (data, SST);
  //if ( Msh()->give_rslts_solver() != SOL_OOFEM ) errol;
  
  if (node && this->give_rslt_NIP(SOL_OOFEM) !=1)  errol;
  
  if (SST == SST_3d) {
    switch (rvtype) {
    case RVTstrn_displ:  check_rslts(step, RTE_global_strain);  data->beCopyOf(this->give_results_dm(step, RTE_global_strain)->give_ptr2val(), 6);  break;
    case RVTstrs_displ:  check_rslts(step, RTE_global_stress);  data->beCopyOf(this->give_results_dm(step, RTE_global_stress)->give_ptr2val(), 6);  break;
    default:             errol;
    }
  }
  else errol;
  
  if (!node)  return 0;
  
  //return this->give_weight();
  _warningg ("volume of brick is approximated");
  
  return this->give_centercoords()->dist2_to(node->give_coords());
}

void Brick :: print_row_VTX (char *str) const
{
  sprintf (str, " %ld %ld %ld %ld %ld %ld %ld %ld",
       ((Node*)points[0])->give_lid_id(0),
       ((Node*)points[1])->give_lid_id(0),
       ((Node*)points[2])->give_lid_id(0),
       ((Node*)points[3])->give_lid_id(0),
       ((Node*)points[4])->give_lid_id(0),
       ((Node*)points[5])->give_lid_id(0),
       ((Node*)points[6])->give_lid_id(0),
       ((Node*)points[7])->give_lid_id(0));
//  sprintf (str, " %ld %ld %ld %ld %ld %ld %ld %ld",
//     ((Node*)points[3])->give_lid_id(0),
//     ((Node*)points[2])->give_lid_id(0),
//     ((Node*)points[1])->give_lid_id(0),
//     ((Node*)points[0])->give_lid_id(0),
//     ((Node*)points[7])->give_lid_id(0),
//     ((Node*)points[6])->give_lid_id(0),
//     ((Node*)points[5])->give_lid_id(0),
//     ((Node*)points[4])->give_lid_id(0));
}


} // namespace midaspace