point.cpp

Go to the documentation of this file.

#include "point.h"

#include "attribute.h"

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

#include <math.h>


namespace midaspace {

void Point :: attributes_allocation (char attflag, const Attributes *masterat)
{
  if (attributes)  errol;
  if (attflag == '-')  errol; //attributes = NULL;
  
  if (masterat) {
    if      (attflag == 'h')  attributes = new   HNAttribs(this, (PointAttribs*)masterat);
    else if (attflag == 'r')  attributes = new  RANAttribs(this, (PointAttribs*)masterat);
    else if (attflag == 'n')  attributes = new PointAttribs(this,               masterat);
    else _errorr0;
  }
  else {
    PointAttribs *na = NULL;
    long prop = this->give_property();
    if (prop)  na = Pd->give_NdAt_with_prop(prop);
    
    if (na && na->give_classid() == classHNAttribs)  { if (attflag != 'n') errol;  else attflag = 'h'; }
    if (na && na->give_classid() == classRANAttribs) { if (attflag != 'n') errol;  else attflag = 'r'; }
    
    if      (attflag == 'h')  attributes = new   HNAttribs(this, (long)0, NULL, na);
    else if (attflag == 'r')  attributes = new  RANAttribs(this, (long)0, NULL, na);
    else if (attflag == 'n')  attributes = new PointAttribs(this, (long)0, NULL, na);
    else _errorr0;
  }
}
DOFsPerNode Point :: give_DOFspnod (void) const
{
  if (this->pointAttribs()->give_DOFbc() == NULL || this->pointAttribs()->give_DOFbc()->give_mask() == NULL)
    return Pd->give_DOFspnod();
  
  return DPN_Void;
}

Point :: Point (const Geometry *owner, long gid, const PoinT *coo, char attflag) : GeometryComponent (owner, gid, gid, 0)
{
  if (coo==NULL)  coords.zero();
  else            coords.copy(coo);
  
  superedges.initialize(0,2);
  superfaces.initialize(0,2);
  superelems.initialize(0,2);
  
  duplmaster = NULL;
  
  subdom = -1;
  
  if (attflag != '-')  this->attributes_allocation(attflag, NULL);
}
Point :: Point (const Point *src, char attflag) : GeometryComponent (src)
{
  coords.copy(&src->coords);
  subdom    = src->subdom;
  
  superedges.be_copy_of(&src->superedges);
  superfaces.be_copy_of(&src->superfaces);
  superelems.be_copy_of(&src->superelems);
  
  duplmaster = NULL;  if (src->duplmaster) errol;
  
  if (attflag != '-') {
    const Attributes *a = (const Attributes*)src->give_pointAttribs();  if (!a) errol;
    this->attributes_allocation(attflag, a);
  }
}
Point :: ~Point()
{
}

void Point :: initialize (void)
{
  GeometryComponent :: initialize();
}
void Point :: finitialize (void)
{
  GeometryComponent :: finitialize();
  
  this->mg_connectivity_assembling();
  
  long i;
  long hinge = pointAttribs()->give_hinge();
  pointAttribs()->set_hinge(0);
  
  switch (hinge) {
  case 1:  // full hinge
    //for (i=0; i<superelems(); i++)
    //  ((Element*)superelems[i])->set_fullhinge_at(this);
    errol;
    break;
  case 2: { // cross hinge
    if (superelems() != 4)  errol;
    Beam *b1 = (Beam*)dynamic_cast<const Beam*>(superelems[0]);  if (!b1) errol; // if (b11->give_nno() != 2) errol;
    Beam *b2 = (Beam*)dynamic_cast<const Beam*>(superelems[1]);  if (!b2) errol; // if (b12->give_nno() != 2) errol;
    Beam *b3 = (Beam*)dynamic_cast<const Beam*>(superelems[2]);  if (!b3) errol; // if (b21->give_nno() != 2) errol;
    Beam *b4 = (Beam*)dynamic_cast<const Beam*>(superelems[3]);  if (!b4) errol; // if (b22->give_nno() != 2) errol;
    
    if (false) {
      VectoR v1, v2, v3, v4;
      b1->give_vector(&v1);
      b2->give_vector(&v2);
      b3->give_vector(&v3);
      b4->give_vector(&v4);
      
      if (v1.is_parallel_with(&v2, ZERO)) {
    if (!v3.is_parallel_with(&v4, ZERO)) _errorr("non straight beam in cross-hinge");
    if ( v1.is_parallel_with(&v3, ZERO)) _errorr("two straight beam in cross-hinge are parallel");
      }
      else {
    if (v1.is_parallel_with(&v3, ZERO)) {
      if (!v2.is_parallel_with(&v4, ZERO)) _errorr("non straight beam in cross-hinge");
      Beam *aux = b2;
      b2 = b3;
      b3 = aux;
    }
    else {
      if (!v1.is_parallel_with(&v4, ZERO)) _errorr("non straight beam in cross-hinge");
      if (!v2.is_parallel_with(&v3, ZERO)) _errorr("non straight beam in cross-hinge");
      Beam *aux = b2;
      b2 = b4;
      b4 = aux;
    }
      }
    }
    if (true) {
      const CrossSection *cs1 = b1->give_elemAttribs()->give_cs();
      if (cs1 == b2->give_elemAttribs()->give_cs()) {
    if (cs1 == b3->give_elemAttribs()->give_cs())  _errorr("cross-hinge: mismatch in cross-sections");
      }
      else {
    if (cs1 == b3->give_elemAttribs()->give_cs()) {
      Beam *aux = b2;
      b2 = b3;
      b3 = aux;
    }
    else {
      if (cs1 != b4->give_elemAttribs()->give_cs())  _errorr("cross-hinge: mismatch in cross-sections");
      Beam *aux = b2;
      b2 = b4;
      b4 = aux;
    }
      }
      
      if (b3->give_elemAttribs()->give_cs() != b4->give_elemAttribs()->give_cs())  _errorr("cross-hinge: mismatch in cross-sections");
    }
    
    Node *node1 = dynamic_cast<Node*>(this);  if (!node1) errol;
    Node *node2 = new HangingNode(node1);
    
    double aux = 1.0;
    if (this->give_DOFspnod() != DPN_DxyzRxyz) errol;
    int n = 6;
    int val[6];
    val[0] = val[1] = val[2] = 2;
    val[3] = val[4] = val[5] = 0;
    
    ((HNAttribs*)node2->pointAttribs())->initialize_general(1, (const Node**)&node1, &aux, n, val);
    
    ((Geometry*)Geom)->add_another_Pjnt(node2);
    
    node1->delete_connectivity();
    node2->delete_connectivity();
    
    b3->switch_node_pointer_in_all_components(node1, node2, false);
    b4->switch_node_pointer_in_all_components(node1, node2, false);
    
    b1->connectivity_assembling(true);
    b2->connectivity_assembling(true);
    b3->connectivity_assembling(true);
    b4->connectivity_assembling(true);
    
    break;
  }
  default: if (hinge > 0) errol;
  }
  
  DOFsPerNode dpnfn = this->give_DOFspnod();    // dofs per node from node
  
  if (dpnfn == DPN_Void) {
    for (i=0; i<superelems(); i++)
      if (superelems[i]->give_elemAttribs()->give_conDOFs(this))
    errol;
  }
  else {
    // breakne se, kdyz element ma stejny dpnfn && nema conDOF => pak podpira plne uzel a vse je v poradku
    for (i=0; i<superelems(); i++)
      if (superelems[i]->give_elemAttribs()->give_conDOFs(this) == NULL  &&  superelems[i]->give_elemAttribs()->give_dpn() == dpnfn)
    break;
    
    // nebylo breaknuto, tak to je hinge
    if (i == superelems()) {
      DOFsPerNode dpne;    // dofs per node from element
      // kontrola, ze to je skutecne full hinge
      // MOMENTALNE ZAKOMENTOVANY, U 2D PRIKLADUS SE STANE ZE give_conDOFs(this) VRATI NULL, TIM PADEM CHYBA
      for (i=0; i<superelems(); i++) {
    dpne = superelems[i]->give_elemAttribs()->give_dpn();
    if ( dpne != DPN_DxyzR___  &&  ( dpne != DPN_DxyzRxyz  || (! superelems[i]->give_elemAttribs()->give_conDOFs(this)->is_full_hinge()) ) )
      errol;
    
    pointAttribs()->reset_dofbc_support_all_rot_dofs();
      }
    }
  }
  
  
  //pointAttribs()->finitialize();
}
void Point :: checkConsistency (void) const
{
  GeometryComponent :: checkConsistency();
  
  // attributes
  //attributes->checkConsistency();
}

const Edge* Point :: give_superedge (long nn, const Point **nods) const
{
  for (int i=0; i<superedges(); i++)
    if (superedges[i]->give_nno() == nn)
      if ( superedges[i]->give_points()->has_these_members(nn, nods) )
    return superedges[i];
  
  return NULL;
}

long Point :: give_property (void) const
{
  if (mproperty.give_size()) { if (mproperty.give_size() > 1)  errol;  return mproperty[0]; }
  if (eproperty.give_size()) { if (eproperty.give_size() > 1)  errol;  return eproperty[0]; }
  if (fproperty.give_size()) { if (fproperty.give_size() > 1)  errol;  return fproperty[0]; }
  if (vproperty.give_size()) { if (vproperty.give_size() > 1)  errol;  return vproperty[0]; }
  return -1;
}
long Point :: give_property (int dim, long id) const
{
  switch (dim) {
  case 0:  if (! mproperty.give_size() > id)  errol;  return mproperty[id];  break;
  case 1:  if (! eproperty.give_size() > id)  errol;  return eproperty[id];  break;
  case 2:  if (! fproperty.give_size() > id)  errol;  return fproperty[id];  break;
  case 3:  if (! vproperty.give_size() > id)  errol;  return vproperty[id];  break;
  default: _errorr("unsupported dimension");
  }
  return -1;
}
void Point :: set_property (int dim, long val)
{
  switch (dim) {
  case 0:  if (!mproperty.give_size()) { mproperty.resize_ignore_vals(1);  mproperty[0] = val; }  else _errorr("main property set already");  break;
  case 1:  if (!eproperty.give_size()) { eproperty.resize_ignore_vals(1);  eproperty[0] = val; }  else _errorr("edge property set already");  break;
  case 2:  if (!fproperty.give_size()) { fproperty.resize_ignore_vals(1);  fproperty[0] = val; }  else _errorr("face property set already");  break;
  case 3:  if (!vproperty.give_size()) { vproperty.resize_ignore_vals(1);  vproperty[0] = val; }  else _errorr("volm property set already");  break;
  default: _errorr("unsupported dimension");
  }
}
void Point :: reset_property (int dim, long val)
{
  switch (dim) {
  case 0:  mproperty.resize_ignore_vals(1);  mproperty[0] = val;  break;
  case 1:  eproperty.resize_ignore_vals(1);  eproperty[0] = val;  break;
  case 2:  fproperty.resize_ignore_vals(1);  fproperty[0] = val;  break;
  case 3:  vproperty.resize_ignore_vals(1);  vproperty[0] = val;  break;
  default: _errorr("unsupported dimension");
  }
}
void Point :: add_property (int dim, long val)
{
  switch (dim) {
  case 0:  mproperty.add(val);  break;
  case 1:  eproperty.add(val);  break;
  case 2:  fproperty.add(val);  break;
  case 3:  vproperty.add(val);  break;
  default: _errorr("unsupported dimension");
  }
}



// ********************************************************
// ***   DUPLICITY   **************************************
// ********************************************************
//
void Point :: make_invisible (Point *master, bool duplcheck)
{
  check_connectivity();
  
  int i;
  for (i=0; i<superedges(); i++)  const_cast<Edge*   >(superedges[i])->switch_node_pointer (this, master, duplcheck);
  for (i=0; i<superfaces(); i++)  const_cast<Face*   >(superfaces[i])->switch_node_pointer (this, master, duplcheck);
  for (i=0; i<superelems(); i++) {
    Element *elem = const_cast<Element*>(superelems[i]);
    elem->switch_node_pointer (this, master, duplcheck);
  }
  
  ((Problem*)Pd)->switch_node_pointer (this, master);
}

void Point :: integrate_duplicated_point (Point *slave)
{
  if (this->duplmaster)  _errorr ("this is slave already");
  
  // zatim tady davam ze to musi byt equal, pac nevim jestli mam akumulovat attributy nebo ne
  this->give_pointAttribs()->checkEqualityWith(slave->give_pointAttribs());
  
  superedges.add_unique (slave->give_superedges()[0]);
  superfaces.add_unique (slave->give_superfaces()[0]);
  superelems.add_unique (slave->give_superelems()[0]);
  //udelej ascendent dodelat
}

bool Point :: invisible_duplicated (char flag)
{
  // flag is here for compatibility with Cell::invisible_duplicated
  
  // no duplicity
  if (duplmaster == NULL)  return false;
  
  // duplicated
  this->duplmaster->integrate_duplicated_point(this);
  this->make_invisible(duplmaster, true);
  
  return true;
}


void Point :: print_row_VTK (FILE *stream) const
{
  // local coords
  fprintf (stream, "% .6e % .6e % .6e\n", coords[0], coords[1], coords[2]);
}
void Point :: print_row_VTX (char *str) const
{
  // local coords
  sprintf (str, "% .6e % .6e % .6e", coords[0], coords[1], coords[2]);
}


Vertex :: Vertex (const Geometry *owner, long gid, const PoinT *coo, char attflag) : Point (owner, gid, coo, attflag)
{
}
//)

Vertex :: ~Vertex()
{
  
}

void Vertex :: initialize (void)
{
  Point :: initialize();
  
  long hinge = pointAttribs()->give_hinge();
  if (hinge == 1) {
    this->mg_connectivity_assembling();
    
    pointAttribs()->set_hinge(0);
    
    for (long i=0; i<superelems(); i++)
      ((Element*)superelems[i])->set_fullhinge_at(this);
  }
  
}

void Vertex :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  switch (fff) {
  case FFF_T3d:
    // name + id
    fprintf (stream, "vertex %3ld   xyz ", this->ID+1);
    
    // local coords
    for (int i=0; i<3; i++)  fprintf (stream, " %21.14e", this->give_coord(i));
    
    // property
    fprintf (stream, "   property %ld", this->ID+1);
    break;
  default: errol;
  }
  
  // end
  if (endline)  fprintf (stream, "\n");
}


Node :: Node (const Geometry *owner, long gid, const PoinT *coo, char attflag) : Point (owner, gid, coo, attflag)
{
  cdom = 0;  lid = NULL;
  mdl_masterex = NULL;
  mdl_masterel = NULL;
  mdl_masteref = NULL;
  resultsN = NULL;
  resultsE = NULL;
}

Node :: Node (const Vertex *src, char attflag) : Point (src, attflag)
{
  cdom = 0;           //if (src->cdom   ) errol;
  lid = NULL;         //if (src->lid    ) errol;
  mdl_masterex = NULL;
  mdl_masterel = NULL;
  mdl_masteref = NULL;
  resultsN = NULL;    //if (src->results) errol;
  resultsE = NULL;    //if (mdl_mstr->results) errol;
}
Node :: Node (const Node *src, char attflag) : Point (src, attflag)
{
  cdom = 0;            if (src->cdom        ) errol;
  lid = NULL;          if (src->lid         ) errol;
  mdl_masterex = NULL; //if (src->mdl_masterex) errol;
  mdl_masterel = NULL; //if (src->mdl_masterel) errol;
  mdl_masteref = NULL;
  resultsN = NULL;     if (src->resultsN    ) errol;
  resultsE = NULL;     if (src->resultsE    ) errol;
}
Node :: ~Node()
{
  delete [] lid;
  deallocateCheckUno (resultsN, Msh()->give_rslts_nsteps(), cRTN, false);
  
  if (resultsE) {
    for (long i=0; i<Msh()->give_nRegions(); i++)
      if (resultsE[i])
    deallocateCheckUno (resultsE[i], Msh()->give_rslts_nsteps(), cRTE, false);
    
    delete [] resultsE;
  }
}

void Node :: initialize (void)
{
  Point :: initialize();
}
void Node :: checkConsistency (void) const
{
  Point :: checkConsistency();
  
  // *** check of position node in domains ***
  if ( cdom == 1  ||  cdom > NumDomains() )
    _errorr ("Mishmash in domains");
  
# ifdef DEBUG
# endif
}

const LCS_pure* Node :: give_mdl_lcs (void)
{
  if (mdl_masterel) {
    if (mdl_masterel->give_dimension() == 2)
      return ((Gelement*)mdl_masterel)->give_cg()->give_lcs();
    else if (mdl_masterel->give_dimension() == 1) {
      const Edge* edge = mdl_masterel->give_edge(0);
      if (edge->give_numsuperface()                    != 1) errol;
      if (edge->give_superface(0)->give_numsuperelem() != 1) errol;
      return ((Element*)edge->give_superface(0)->give_superelem(0))->give_cg()->give_lcs();
    }
  }
  
  if (mdl_masterex) {
    if (mdl_masterex->give_numsuperface()                    != 1) errol;
    if (mdl_masterex->give_superface(0)->give_numsuperelem() != 1) errol;
    
    return ((Element*)mdl_masterex->give_superface(0)->give_superelem(0))->give_cg()->give_lcs();
  }
  
  errol;
  return NULL;
}



void Node :: initialize_domli (bool shared, long dom, long li)
{
  if (!Parallel())  errol;
  
  // cdom, lid
  if (shared){ cdom =    1; lid = new long[NumDomains()]; fill_all_by (lid, NumDomains(), (long)-1); lid[dom] = li; }  // node is     shared
  else       { cdom = -dom; lid = new long[     1    ];                                              lid[ 0 ] = li; }  // node is not shared
}

void Node :: initialize_parallel (long nd, long *dom)
{
# ifdef DEBUG
  if (cdom!=0 || lid!=NULL)  _errorr ("ERROR");
# endif
  
  cdom = nd;
  if      (nd<1  && dom==NULL) { lid = new long[1];  lid[0] = ((Mesh*)Msh())->NumNodes_plus(-nd); }
  else if (nd>1  && dom!=NULL) { lid = new long[NumDomains()];
    fill_all_by (lid, NumDomains(), (long)-1);
    for (int i=0; i<nd; i++) lid[dom[i]] = ((Mesh*)Msh())->NumNodes_plus(dom[i]);
  }
  else _errorr0;
}

void Node :: add_domain (long gid, const char *&str, femFileFormat fff, bool shared, long dom, long li)
{
  switch (fff) {
  case FFF_T3d: {
    PoinT coo;
    coo.scan_xyz (str);
    if (!coords.is_identical_to(&coo))  _errorr ("coords differ");
    
    long parentdim, trash, prop;
    
    SP_scan_number (str, parentdim);
    SP_scan_number (str, trash);
    SP_scan_number (str, prop);
    
    if (prop != give_property (EntityType_type2dim_jktkT3d (parentdim), 0))  _errorr ("property differs");
    
    break;
  }
  default: _errorr("unsupported fem file format");
  }
  
  if (!Parallel())  _errorr ("There no lid for non-Parallel analysis");
  if (!shared)      _errorr ("The node has to be shared");
  if (cdom<=0) {
    // toto odstranuje Danovu neschopnost tisknout u prvniho file glob cisla shareovanych uzlu s minusem 
    if (dom==1 && li<4) {
      gid=lid[0]; delete lid;
      lid = new long[NumDomains()];  fill_all_by (lid, NumDomains(), (long)-1);
      lid[-cdom] = gid; cdom = 1;
    }
    else _errorr ("The node has to be shared, cdom<=0");
  }
  lid[dom] = li;
  cdom++;
}


bool Point :: is_on_3d_element (void) const
{
  for (long i=0; i<superelems(); i++)
    if (superelems[i]->give_dimension() == 3)
      return true;
  
  return false;
}

// prida existujici node (nemusi se allocovat) do domeny == rozsiri se lid[]
void Node :: add_me_to_domain (long did)
{
  if (!Parallel())  _errorr ("There no lid for non-Parallel analysis");
  
  if (this->give_lid(did) != -1)  return;
  
  if (cdom<1) {
    long aux = lid[0];
    delete [] lid;
    lid = new long [NumDomains()];
    for (int i=0;i<NumDomains();i++)  lid[i] = -1;
    lid[-cdom] = aux;
    cdom = 1;
  }
  cdom++;
  lid[did] = ((Mesh*)Msh())->NumNodes_plus(did);
}

// find parent subdomain
bool Node :: find_parent_subdom (long sdid, long *nn, long *level)
{
  if ( subdom >= 0 ) {
    if (level[0]  &&  subdom != sdid)  _errorr4("Error in subdom, %ld != %ld, ID = %ld", subdom, sdid, ID);
    return false;   // subdom flag already set
  }
  
  level[0]++;
  
  if (level[0] > 50000)  return true;  // allowed number of recursive loops
  
  nn[0]++;
  subdom = sdid;
  
  for (long i=0; i<superelems(); i++)
    for (long j=0; j<superelems[i]->give_nno(); j++)
      if ( ((Node*)superelems[i]->give_point(j))->find_parent_subdom (sdid,nn,level) ) {
    nn[0]--;
    subdom = -2;
    level[0]--;
    return true;
      }
  
  level[0]--;
  return false;
}

//double Node :: give_absDisp (void) const
//{
//  return sqrt (displ[0]*displ[0] + displ[1]*displ[1] + displ[2]*displ[2]);
//}


double Node :: give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step)
{
  double w, elw;
  Dvctr eldata;
  
  w = 0.0;
  data->resize_ignore_vals(6);
  data->zero();
  for (long i=0; i<superelems(); i++) {
    elw = ((FElement*)superelems[i])->give_ssstate (&eldata, SST, rvtype, type, step, this);
    if (eldata.give_size() == 6) {
      w += elw;
      data->addtms(eldata, elw);
    }
  }
  if(w>1.0e-10)
    data->dvdby(w);
  return w;
}


void Node :: set_master_component (long prop, const Model *model, int parenttype)
{
  if (parenttype == 1)                 mdl_masterex = dynamic_cast<const Vertex   *>(model->give_Pjnt(prop-1));
  else if (prop <= model->give_Nels()) mdl_masterel = dynamic_cast<const Gelement *>(model->give_Elem(prop-1));
  else {
    switch (EntityType_type2dim_jktkT3d (parenttype)) {
    case 1:  mdl_masteref = model->give_Edge(prop - model->give_Nels() - 1);  break;
    case 2:  mdl_masteref = model->give_Face(prop - model->give_Nels() - 1);  break;
    default: errol;
    }
  }
# ifdef DEBUG
  if (mdl_masterex == NULL && mdl_masterel == NULL && mdl_masteref == NULL)  errol;
# endif
}

const PointDOFsBCPM* give_uniq_dofbc (int dim, const GPA<const Element> *superelems, const char *name, long id)
{
  const PointDOFsBCPM *dofbc, *dofbc2;
  
  int i;
  long supid;
  dofbc = NULL;
  
  for (; dim<=3; dim++) {
    for (i=0; i<(*superelems)(); i++) {
      if ( dim != (*superelems)[i]->give_dimension() )  continue;
      
      dofbc2 = ((const Gelement*)(*superelems)[i])->give_gelemAttribs()->give_dofbc();
      if (dofbc2) {
    if (dofbc==NULL) { dofbc = dofbc2;  supid = (*superelems)[i]->give_ID(); }
    else
      if (!dofbc->is_identical_with(dofbc2))
        _errorr5("Different boundary conditions at elements %ld and %ld meet at %s %ld", supid+1, (*superelems)[i]->give_ID()+1, name, id+1);
      }
    }
    
    if (dofbc) break;
  }
  
  return dofbc;
}

void Node :: read_input (const char *&str, femFileFormat fff)
{
  switch (fff) {
  case FFF_T3d: {
    coords.scan_xyz (str);
    
    long parenttype, trash, prop;
    
    SP_scan_number (str, parenttype);
    SP_scan_number (str, trash);
    SP_scan_number (str, prop);
    
    if (prop)  this->set_property(EntityType_type2dim_jktkT3d (parenttype), prop);
    
    // master model element
    const Model *model = ((Mesh*)Geom)->give_Master_model();
    
    // master model element
    // toto musim predelat na kontrolu vsech property, az mi dan zacne davat oproperty
    // 
    if (model) {
      this->set_master_component (prop, model, parenttype);
      this->attributes_allocation ('n', mdl_masterex ? mdl_masterex->give_attributes() : (mdl_masterel ? mdl_masterel->give_attributes() : NULL));
    }
    else
      this->attributes_allocation('n', NULL);
    
    //* inherit boundary condition
    if (model) {
      const PointDOFsBCPM *dofbc = this->pointAttribs()->give_DOFbc();
      
      if (dofbc) break;
      
      if (mdl_masterex)
    dofbc = give_uniq_dofbc (1, mdl_masterex->give_superelems(), "vertex", mdl_masterex->give_ID());
      else if (mdl_masterel) {
    const Facedge* fe = mdl_masterel->give_same_dimension_facedge();
    dofbc = give_uniq_dofbc (fe->give_dimension(), fe->give_superelems(), "element", mdl_masterel->give_ID());
      }
      else
    dofbc = give_uniq_dofbc (mdl_masteref->give_dimension(), mdl_masteref->give_superelems(), "edge/face", -1);
      
      if (dofbc)
    this->pointAttribs()->set_dofbc_copy_of(dofbc);
    }
    break;
  }
  case FFF_OOFEM:
    SP_scan_expected_word_exit   (str, "coords", "Invalid structure of the given OOFEM input file", CASE);
    SP_scan_expected_number_exit (str,        3, "Invalid structure of the given OOFEM input file, the number of coords");
    
    coords.scan_xyz (str);
    
    pointAttribs()->initialize_from (str, fff);
    break;
  case FFF_SIFEL:
    coords.scan_xyz (str);
    
    pointAttribs()->initialize_from (str, fff);
    break;
  case FFF_JKTK:
    coords.scan_xyz (str);
    
    long i, auxl, et, prop;
    SP_scan_number (str, auxl);
    for (i=0; i<auxl; i++) {
      SP_scan_number (str, et);
      SP_scan_number (str, prop);
      if (prop) this->add_property (EntityType_type2dim_jktkT3d (et), prop);
    }
    break;
  default: _errorr("unsupported fem file format");
  }
}


void Node :: read_output_OOFEM (FILE *stream, long step, ResultTypesAtNode rt)
{
  const PointDOFsBCPM* bc = pointAttribs()->give_DOFbc();
  if (rt == RTN_reaction && (bc == NULL  ||  bc->is_supported() == false))  return;
  
  const char *str;
  char LINE[1023];
  
  int ndofs = this->give_pointAttribs()->give_nDOFs();
  Dvctr *values = new Dvctr(ndofs);
  
  if (rt == RTN_displacement) {
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word3_exit  (str, "Node", "RigidArmNode", "HangingNode", "Invalid structure of OOFEM output file", CASE);
    SP_scan_expected_number_exit (str, ID+1, "Invalid structure of OOFEM output file");
    
    //
    for (int i=0; i<ndofs; i++) {
      fgets (LINE, 1023, stream);  str=LINE;
      // dof
      SP_scan_expected_word_exit   (str, "dof", "Invalid structure of OOFEM output file", CASE);
      SP_scan_expected_number_exit (str,   i+1, "Invalid structure of OOFEM output file");
      // d
      SP_scan_expected_word_exit   (str,   "d", "Invalid structure of OOFEM output file", CASE);
      SP_scan_number (str, (*values)[i]);
    }
  }
  else if (rt == RTN_reaction) {
    values->zero();
    
    for (int i=0; i<ndofs; i++) {
      if (bc->give_att(i)) {
      fgets (LINE, 1023, stream);  str=LINE;
      SP_scan_expected_word_exit   (str, "Node"    , "Invalid structure of OOFEM output file", CASE);
      SP_scan_expected_number_exit (str, this->ID+1, "Invalid structure of OOFEM output file");
      SP_scan_expected_word_exit   (str, "iDof"    , "Invalid structure of OOFEM output file", CASE);
      SP_scan_expected_number_exit (str,        i+1, "Invalid structure of OOFEM output file");
      SP_scan_expected_word_exit   (str, "reaction", "Invalid structure of OOFEM output file", CASE);
      
      SP_scan_number               (str, (*values)[i]);
      
      SP_scan_expected_word_exit   (str, "[bc-id:" , "Invalid structure of OOFEM output file", CASE);
      SP_scan_expected_word_exit   (str, "1]"      , "Invalid structure of OOFEM output file", CASE);
      }
    }
  }
  else errol;
  
  this->add_resultN (values, step, rt);
}

void Node :: read_output_SIFEL (FILE *stream, long step, ResultTypesAtNode rt)
{
  const PointDOFsBCPM* bc = pointAttribs()->give_DOFbc();
  if (rt == RTN_reaction && (bc == NULL  ||  bc->is_supported() == false))  return;
  
  const char *str;
  char LINE[1023];
  
  if (rt == RTN_displacement) {
    int ndofs = this->give_pointAttribs()->give_nDOFs();
    
    Dvctr *values = new Dvctr(ndofs);
    
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit   (str, "Node", "Invalid structure of SIFEL output file", CASE);
    SP_scan_expected_number_exit (str, ID+1,   "Invalid structure of SIFEL output file");
    
    //
    char r[6];
    for (int i=0; i<ndofs; i++) {
      sprintf(r,"r_%d=",i+1);
      SP_scan_expected_word_exit   (str, r, "Invalid structure of SIFEL output file", CASE);
      SP_scan_number (str, (*values)[i]);
    }
    this->add_resultN (values, step, rt);
  }
  else if (rt == RTN_transpunknowns) {
    int nmed = this->give_pointAttribs()->give_nDOFs();
    Dvctr *values = new Dvctr(nmed);
    
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit   (str, "Node", "Invalid structure of SIFEL output file", CASE);
    SP_scan_expected_number_exit (str, ID+1,   "Invalid structure of SIFEL output file");
    
    //
    char r[6];
    for (int i=0; i<nmed; i++) {
      sprintf(r,"r_%d=",i+1);
      SP_scan_expected_word_exit   (str, r, "Invalid structure of SIFEL output file", CASE);
      SP_scan_number (str, (*values)[i]);
    }
    this->add_resultN (values, step, rt);
  }
  else if (rt == RTN_reaction) {
    int ndofs = this->give_pointAttribs()->give_nDOFs();
    
    Dvctr *values = new Dvctr(ndofs);
    values->zero();
    
    fgets (LINE, 1023, stream);  str=LINE;
    SP_scan_expected_word_exit   (str, "Node", "Invalid structure of SIFEL output file", CASE);
    SP_scan_expected_number_exit (str, ID+1,   "Invalid structure of SIFEL output file");
    
    //
    char r[6];
    for (int i=0; i<ndofs; i++) {
      sprintf(r,"R_%d=",i+1);
      SP_scan_expected_word_exit   (str, r, "Invalid structure of SIFEL output file", CASE);
      SP_scan_number (str, (*values)[i]);
      if (!bc->give_att(i)) {
    if ((*values)[i] > 1.e-12)  errol;
    else                        (*values)[i] = 0.0;
      }
    }
    this->add_resultN (values, step, rt);
  }
  else errol;
}

void Node :: print_row (FILE *stream, femFileFormat fff, bool endline, long did) const
{
  char flag = 'a'; // to je tady historicky, posilalo se to jako argument, ale ted to nikdo nevola, tak jsem to zrusil
  long i;
  
  if (fff == FFF_OOFEM) {
    if (flag == 'a' || flag == 's') {
      // name + id
      fprintf (stream, "%s %4ld", this->give_node_name(), this->give_lid_id(did)+1);
      
      // local coords
      fprintf (stream, "  coords 3");
      for (i=0; i<3; i++)  fprintf (stream, " %16.9e", this->give_coord(i));
      
      // global id
      if (Parallel())  fprintf (stream, "  globNum %4ld", this->give_ID()+1);
      
      // boundary conditions
      this->pointAttribs()->print_row(stream, fff, did);
    }
    if (flag == 'a' || flag == 'e') {
      // parallel, shared
      if (this->is_shared()) {
    fprintf (stream, "  Shared partitions %4ld", this->give_cdom()-1);
    for (i=0; i<NumDomains(); i++)
      if (this->give_lid(i)!=-1 && i!=did)
        fprintf (stream," %4ld",i);
      }
    }
  }
  else if (fff == FFF_SIFEL) {
    fprintf (stream, "%6ld", this->give_ID()+1);
    
    // local coords
    for (i=0; i<3; i++)  fprintf (stream, " %13.6e", this->give_coord(i));
    
    // boundary conditions
    this->pointAttribs()->print_row(stream, fff, did);
  }
  else if (fff == FFF_ANSYS) {
    fprintf (stream, "N,%6ld", this->give_ID()+1);
    
    // local coords
    for (i=0; i<3; i++)  fprintf (stream, ", %13.6e", this->give_coord(i));
  }
  else if (fff == FFF_JKTK) {
    if (Pd->give_PDBO(PDBO_melnik)) {
      fprintf (stream, "%ld", this->give_ID()+1);
      // local coords
      for (i=0; i<3; i++)  fprintf (stream, " %16.10e", this->give_coord(i));
      // properties
      long nprop = mproperty.give_size() + eproperty.give_size() + fproperty.give_size() + vproperty.give_size();
      //if (nprop == 0) _errorr2("There is no property at node %ld, required by JKTK file format", this->give_ID()+1);
      fprintf (stream, " %ld", nprop);
      for (i=0; i<vproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(3), vproperty[i]);
      for (i=0; i<fproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(2), fproperty[i]);
      for (i=0; i<eproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(1), eproperty[i]);
      for (i=0; i<mproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(0), mproperty[i]);
      
      // if (this->give_mproperty_cnt())  fprintf (stream, "  1  4 %ld", this->give_propertyORzero());
      // else if (oproperty.give_size()) {
      //    fprintf (stream, " %ld", oproperty.give_size()/2);
      //    for (i=0; i<oproperty.give_size(); i+=2)
      //      fprintf (stream, " %ld %ld", oproperty[i], oproperty[i+1]);
      // }
    }
    else {
      fprintf (stream, "%8ld", this->give_ID()+1);
      // local coords
      for (i=0; i<3; i++)  fprintf (stream, " %13.6e", this->give_coord(i));
      // properties
      // properties
      long nprop = (mproperty.give_size() ? mproperty.give_size() : 1) + eproperty.give_size() + fproperty.give_size() + vproperty.give_size();
      if (nprop == 0) _errorr2("There is no property at node %ld, required by JKTK file format", this->give_ID()+1);
      fprintf (stream, " %ld", nprop);
      
      if (mproperty.give_size() == 0)
    fprintf (stream, " %d %d", EntityType_dim2type_jktkT3d(0), 0);
      else
    for (i=0; i<mproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(0), mproperty[i]);
      
      for (i=0; i<eproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(1), eproperty[i]);
      for (i=0; i<fproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(2), fproperty[i]);
      for (i=0; i<vproperty.give_size(); i++)  fprintf (stream, " %d %ld", EntityType_dim2type_jktkT3d(3), vproperty[i]);
    }
  }
  else errol;
  
  fprintf (stream, "\n");
}

void Node :: allocate_resultsN (void)
{
  if (resultsN) return;
  
  if (Msh()->give_rslts_nsteps() == 0)  _errorr("allocate results with nsteps == 0, mayby replace -IN_VTKaddata by -IN_VTKaddataRSLT"); 
  
  resultsN = new Array1d**[Msh()->give_rslts_nsteps()];
  
  for (long i=0; i<Msh()->give_rslts_nsteps(); i++) {
    resultsN[i] = new Array1d*[cRTN];
    memset (resultsN[i], 0, cRTN*sizeof(Array1d*));
  }
}

void Node :: allocate_resultsE (long regid)
{
  if (resultsE && resultsE[regid]) return;
  
  if (resultsE == NULL) {
    resultsE = new Array***[Msh()->give_nRegions()];
    memset (resultsE, 0, Msh()->give_nRegions()*sizeof(Array***));
  }
  
  long nsteps = Msh()->give_rslts_nsteps();  if (nsteps == 0)  _errorr("allocate results with nsteps == 0, mayby replace -IN_VTKaddata by -IN_VTKaddataRSLT"); 
  
  resultsE[regid] = new Array**[nsteps];
  
  while(nsteps--) {
    resultsE[regid][nsteps] = new Array*[cRTE];
    memset (resultsE[regid][nsteps], 0, cRTE*sizeof(Array*));
  }
}

void Node :: add_resultN (               Xvctr *rslt,             long step, ResultTypesAtNode rt) { this->allocate_resultsN();      if (this->resultsN       [step][rt]) errol;  this->resultsN       [step][rt] = rslt; }
void Node :: set_resultN (long s, const double *rslt,             long step, ResultTypesAtNode rt) { this->allocate_resultsN();      if (this->resultsN       [step][rt]) errol;  this->resultsN       [step][rt] = new Dvctr(s, rslt); }
void Node :: set_resultN (        const double  rslt,             long step, ResultTypesAtNode rt) { this->allocate_resultsN();      if (this->resultsN       [step][rt]) errol;  this->resultsN       [step][rt] = new Dscal(   rslt); }
void Node :: add_resultE (               Dvctr *rslt, long regid, long step, ResultTypesAtElem rt) { this->allocate_resultsE(regid); if (this->resultsE[regid][step][rt]) errol;  this->resultsE[regid][step][rt] = rslt; }
void Node :: set_resultE (        const  Dmtrx *rslt, long regid, long step, ResultTypesAtElem rt) { this->allocate_resultsE(regid); if (this->resultsE[regid][step][rt]) errol;  this->resultsE[regid][step][rt] = new Dmtrx(   rslt); }
void Node :: set_resultE (        const VectoR *rslt, long regid, long step, ResultTypesAtElem rt) { this->allocate_resultsE(regid); if (this->resultsE[regid][step][rt]) errol;  this->resultsE[regid][step][rt] = new Dvctr(   rslt); }

const Dscal* Node :: give_resultsN_ds (            long step, ResultTypesAtNode rt) const { if (resultsN       [step][rt] == NULL)  return NULL;  Dscal *rslt = dynamic_cast<Dscal*>(resultsN       [step][rt]);  if (!rslt) errol;  return rslt; }
const Dvctr* Node :: give_resultsN_dv (            long step, ResultTypesAtNode rt) const { if (resultsN       [step][rt] == NULL)  return NULL;  Dvctr *rslt = dynamic_cast<Dvctr*>(resultsN       [step][rt]);  if (!rslt) errol;  return rslt; }
//
const Dvctr* Node :: give_resultsE_dv (long regid, long step, ResultTypesAtElem rt) const { if (resultsE == NULL || resultsE[regid] == NULL) return NULL;
                                                                                            if (resultsE[regid][step][rt] == NULL)  return NULL;  Dvctr *rslt = dynamic_cast<Dvctr*>(resultsE[regid][step][rt]);  if (!rslt) errol;  return rslt; }
const Dmtrx* Node :: give_resultsE_dm (long regid, long step, ResultTypesAtElem rt) const { if (resultsE == NULL || resultsE[regid] == NULL) return NULL;
                                                                                            if (resultsE[regid][step][rt] == NULL)  return NULL;  Dmtrx *rslt = dynamic_cast<Dmtrx*>(resultsE[regid][step][rt]);  if (!rslt) errol;  return rslt; }

void Node :: setup_full_alloc_DOFvals_at (Dvctr *answer, ResultTypesAtNode rt, long step) const
{
  if (rt != RTN_displacement  &&  rt != RTN_reaction)  errol;
  
  const Dvctr *values = this->give_resultsN_dv(step, rt);
  
  if (values) pointAttribs()->setup_full_alloc_DOFvals(answer, values); 
  else        answer->resize_ignore_vals(0);
}

void Node :: give_displcmnt_in_all_steps(Dvctr *data, int indx) const
{
  if (! this->give_resultsN_dv(0, RTN_displacement))  errol;
  
  for (long i=0; i<Msh()->give_rslts_nsteps(); i++)
    data[0][i] = (*(Dvctr*)this->resultsN[i][RTN_displacement])[indx];
  
}



void HangingNode :: initialize (void)
{
  Node :: initialize();
  
  this->initialize_masters_at_domains();
}
void HangingNode :: checkConsistency (void) const
{
  Node :: checkConsistency();
  
  // *** check attributes ***
  if ( dynamic_cast<HNAttribs *>(attributes) == NULL )  _errorr0;
  
  //
  if (Parallel()) {
    for (int did=0; did<NumDomains(); did++) {
      if (cdom<1)  did = -cdom;
      else         if (lid[did] == -1) continue;
      
      this->give_HNattrb()->check_masters_at_domain(did);
      
      if (cdom<1) break;
    }
  }
}

void HangingNode :: initialize_masters_at_domains (void)
{
  if (!Parallel())  return;
  
  for (int did=0; did<NumDomains(); did++) {
    if (cdom<1)  did = -cdom;
    else         if (lid[did] == -1) continue;
    
    ((HNAttribs*)attributes)->add_masters_to_domain(did);
    
    if (cdom<1) break;
  }
}
void HangingNode :: initialize_hn (const PoinT *cartcoord, const Cell *mc, long dimMC, long ordMC, long countMN, const Node **nodes, const PoinT *natcoord)
{
  coords.copy(cartcoord);
  
  ((HNAttribs*)attributes)->initialize_atelemet (mc, dimMC, ordMC, countMN, nodes, natcoord);
}

//void HangingNode :: read_input_OOFEM (const char *str){}

void HangingNode :: print_one_row_HN (FILE *stream, long hid) const
{
  ((HNAttribs*)attributes)->print_one_row (stream, hid);
  
  if (Parallel()){
    fprintf (stream," gid %ld cdom %3ld lid", ID+Msh()->give_Nn(), cdom); 
    if (cdom<1) fprintf (stream," %ld",lid[0]);
    else for (long i=0; i<NumDomains(); i++) fprintf (stream," %ld", lid[i]);
  }
}


void RigidArmNode :: initialize (void)
{
  Node :: initialize();
  
  //if (property2 != -1) {
  //  _errorr0;
  //  //prpos2 = Pd->give_RAN_prpos (property2);
  //  //if (prpos2 == -1) _errorr3 ("Error give_RAN_prpos, property1 = %d, property2 = %d", property, property2);
  //  //
  //  //if (Pd->give_RAN_prpos(this->property) != -1)
  //  //  Pd->give_NdAt(this->property)->checkEqualityWith (Pd->give_NdAt(this->property2));
  //}
  
  // *** finding shared rigid arm node masters ***
  if (Parallel())
    for (int did=0; did<NumDomains(); did++) {
      if (cdom<1)  did = -cdom;
      else         if (lid[did] == -1) continue;
      
      ((RANAttribs*)attributes)->add_master_to_domain(did);
      
      if (cdom<1) break;
    }
}
void RigidArmNode :: checkConsistency (void) const
{
  Node :: checkConsistency();
  
  // *** check attributes ***
  if ( dynamic_cast<RANAttribs *>(attributes) == NULL )  _errorr0;
  
# ifdef DEBUG
# endif
}

//const Node* RigidArmNode :: give_RANmaster (void) const
//{
//  return ((RANAttribs*)attributes)->give_master();
//}

//void RigidArmNode :: print_row_OOFEM (FILE *stream, long did, char flag) const
//{
//  if (flag != 'a')  _errorr ("error");
//  
//  Node :: print_row_OOFEM(stream, did, 's');
//  
//  
//  //long i, ndofs;
//  // master
//  if (this->property2 == -1)
//    ;
//    //pointAttribs()->print_row_OOFEM (stream, did);
//  else {
//    _errorr0;
//    //const RANode *ran1 = Pd->give_RAN(this->prpos );
//    //const RANode *ran2 = Pd->give_RAN(this->prpos2);
//    //
//    //int m1 = ran1->give_master()->give_lid (did)+1;
//    //int m2 = ran2->give_master()->give_lid (did)+1;
//    //ndofs = Pd->give_NdAt(this->give_property())->give_nDOFs();
//    //int *mstr = new int [ndofs];
//    //int *mmsk = new int [ndofs];
//    //int *dftp = new int [ndofs];
//    //
//    //int aux1,aux2;
//    //for (i=0; i<ndofs; i++) {
//    //  aux1 = ran1->give_mask(i);
//    //  aux2 = ran2->give_mask(i);
//    //  if (aux1*aux2)  _errorr ("RAN mismatch");
//    //  if (aux1) mstr[i] = m1;   else if (aux2) mstr[i] = m2;   else mstr[i] = 0;
//    //  if (aux1) mmsk[i] = aux1; else if (aux2) mmsk[i] = aux2; else mmsk[i] = 0;
//    //  dftp[i] = (aux1+aux2 ? 2 : 0);
//    //}
//    //fprintf (stream,"  masters %ld",    ndofs);  for (i=0; i<ndofs; i++) fprintf (stream," %d", mstr[i]);
//    //fprintf (stream,"  dofType %ld",    ndofs);  for (i=0; i<ndofs; i++) fprintf (stream," %d", dftp[i]);
//    //fprintf (stream,"  masterMask %ld", ndofs);  for (i=0; i<ndofs; i++) fprintf (stream," %d", mmsk[i]);
//    //
//    //delete [] mstr; delete [] mmsk; delete [] dftp;
//  }
//  
//  Node :: print_row_OOFEM(stream, did, 'e');
//}
//RigidArmNode :: RigidArmNode (long gid, const Mesh *mg) : Node (gid, mg)
//{
//  //  prpos2 = -1;
//  //property2 = -1;
//}
//RigidArmNode :: RigidArmNode (long gid, long prop, const Mesh *mg, double *coo, bool shared, long dom, long li) : Node (gid, prop, mg, coo, shared, dom, li)
//{
//  //prpos2 = -1;
//  //
//  if (property>9999999) {
//    _errorr0;
//    //property2 = property%10000;
//    //property  = property/10000;
//  }
//  //else property2 = -1;
//}

} // namespace midaspace