cell.h

Go to the documentation of this file.

#ifndef MIDAS_CELL_H
#define MIDAS_CELL_H

#include "alias.h"
#include "geomcomp.h"
#include "compgeom.h"

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


namespace midaspace {

//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS CELL                          *** *** *** ***
//*  ********************************************************************************************
class Cell : public GeometryComponent //, virtual public ComponentGeometry
{
  // nedokazal jsem to oddedit tak aby to chodilo, tak jsme udelal promennou cg
  //*  *** COMPONENTGEOMETRY ***
 protected:
  ComponentGeometry *cg;
  
 public:
  CellGeometry give_cellGeom           (void)               const { return cg->give_elemGeom();       }  
  int          give_dimension          (void)               const { return cg->give_dimension();      }
  const PoinT* give_centercoords       (void)               const { return cg->give_centercoords();   }
  double       give_circum             (void)               const { return cg->give_circum();         }
  double       give_lav                (void)               const { return cg->give_lav();            }
  double       give_GeomWeight1deg     (void)               const { return cg->give_GeomWeight1deg(); }
  bool         is_point_in_sphere      (const PoinT *point) const { return cg->is_point_in_sphere(point);  }
  double       give_characteristic_size(void)               const { return cg->give_characteristic_size(); }
  
  // docasne kuli lcs
  ComponentGeometry* give_cg           (void)                { return cg;       }
  
  
  
  //*  *** COMPONENTS ***
  GPA<const Point> points;        // array of pointers to points at Cell
  GPA<const Edge>  edges;         // array of pointers to edges  at Cell;  assembled only when this->connectivity_assembled == true
  GPA<const Face>  faces;         // array of pointers to faces  at Cell;  assembled only when this->connectivity_assembled == true
  
  //*  *** CONNECTIVITY ***
  bool connectivity_assembled;    // connectivity stats at cell and its components are assembled
  
  //*  *** DUPLICITY ***
  // duplmaster == -1   ... potentially duplicated
  // duplmaster == this ... collapsed
  // duplmaster != NULL ... duplicated
  // Duplmaster has lower ID then this. The duplicity can be chained.
  // So, the chain is descending! More slaves can have same master.
  Cell *duplmaster; // this Cell has identical nodes with master
  
  
 public:
  Cell (classID mecg, long gid, long oid, long prop, const Geometry *mg, long nn, long ne, long nf);
  Cell (classID mecg, const Cell * src);
  virtual ~Cell();
  
  virtual void checkConsistency (void) const {
    GeometryComponent :: checkConsistency();
    if (duplmaster!=NULL && duplmaster!=((Cell*)-1))  errol;
  }
  void check_connectivity_loc (void) const {
    if (! this->connectivity_assembled)
      _errorr ("Connectivity on cell is not initiated");
  }
  
  
  //*  *** SET ***
  void set_point (long i, const Point *val) { points.assign(i, val); }
  void set_point (long i,       long   pid) { points.assign(i, Geom->give_Pjnt(pid)); }
  void set_edge  (long i, const Edge  *val) { edges.assign (i, val); }
  void set_face  (long i, const Face  *val) { faces.assign (i, val); }
  //*  *** GET ***
  long give_nno (void) const { return points(); }
  long give_ned (void) const { return edges(); }
  long give_nfa (void) const { return faces(); }
  //
  const Point* give_point (long i) const {                           return points[i]; }
  const Edge*  give_edge  (long i) const { check_connectivity_loc(); return edges[i]; }
  const Face*  give_face  (long i) const { check_connectivity_loc(); return faces[i]; }
  //
  const GPA<const Point>* give_points (void) const {                       return &points; } // !!!! zkusit vrati &
  const GPA<const Edge>*  give_edges  (void) const { check_connectivity_loc(); return &edges; }
  const GPA<const Face>*  give_faces  (void) const { check_connectivity_loc(); return &faces; }
  //
  long give_point_indx (const Point* p) const { for (int i=0; i<points(); i++) if (points[i] == p) return i;  errol; return 0; }
  //
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const { errol; return 0; }
  
  //*  *** CONNECTIVITY ***
  virtual void connectivity_assembling (bool re=false) = 0;
  virtual void connectivity_removing (void) = 0;
  virtual void switch_myself_at_connectivity (Cell *master);
  
  
  //protected: proc to nejde ???
  bool switch_edge_pointer (Edge *slave, Edge *master);
  bool switch_face_pointer (Face *slave, Face *master);
  
  
 private:
  void  assure_duplicity_master (void);
  Cell* last_duplicity_master (void) const;
 protected:
  virtual void integrate_duplicated_cell (const Cell *slave);
  virtual void switch_node_pointer (Point *slave, Point *master, bool duplcheck);
 public:
  void setup_duplicity_master (Cell *val);
  bool invisible_duplicated (char flag='a');
  
  bool has_same_geom_with (Cell *slave) const;
  
  bool check_collapse (void) { return ! points.is_unique(); }
  
  bool cross_abscissa_node (const PoinT *a1, const PoinT *a2, long cunn, const Point **unnod, const Point *&nod, double &ksi, PoinT *cp) const;
  int  cross_abscissa_face (const PoinT *a1, const PoinT *a2, long cunf, const Face  **unfac, const Cell *&comp, PoinT *nc,  PoinT *cp);
};



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS FACEDGE                         *** *** *** ***
//*  ********************************************************************************************
class Facedge : public Cell
{
  // empty, only connect face and edge for future features
  // classy zde budou primitivni, napr ord je zatim pouze 1
  
 protected:
  GPA<const Element> superelems;     // array of pointers to superior elements; this array should be sorted!!!
  
  int masterel_dim;
  GPA<const Element> masterels;      // array of pointers to master elements - superior elements with the most low dimension
  
 public:
  Facedge (long gid, long oid, long prop, const Geometry *mg, long nn, long ne, long nf) : Cell (classVoid, gid, oid, prop, mg, nn, ne, nf)
  {
    superelems.initialize(0,2);
    masterel_dim = 7;
  }
  //Facedge () : Cell () {    superfaces.initialize(0,2); }
  virtual ~Facedge() { }
  
  //virtual classID give_classid() const = 0;
  
  
  //*  *** SET ***
  void setadd_superelem (const Element *comp, bool uniquecheck) { if (uniquecheck) superelems.add_unique(comp); else superelems.add(comp);  this->setadd_masterel(comp); }
  void remove_superelem (const Element *comp)                   {                  superelems.free_ptr(comp);  superelems.shake_down(); }
  void setadd_masterel  (const Element *comp);
  
  //*  *** GET ***
  long                   give_numsuperelem (void)   const { check_connectivity(); return  superelems();  }
  const Element*            give_superelem (long i) const { check_connectivity(); return  superelems[i]; }
  const GPA<const Element>* give_superelems (void)  const { check_connectivity(); return &superelems;    }
  
  
  const Element*            give_masterel   (void)  const { check_connectivity(); return masterels[0];  }
  const Element*            give_main_masterel_uniq (void) const;
  
        FacedgeAttribs *give_facedgeAttribs (void)       { return (FacedgeAttribs*)attributes; }
  const FacedgeAttribs *give_facedgeAttribs (void) const { return (FacedgeAttribs*)attributes; }
  

 protected:
  virtual void integrate_duplicated_cell (const Facedge *slave);
  
 public:
  virtual long give_ord (void) const = 0;
  
  
  const GPA<const BoundaryCond>* give_loads  (const Element *fe) const
  {
    return (this->give_masterel() == fe) ? this->give_attributes()->give_loads() : NULL;
  }
  
  long                           give_cloads (const Element *fe) const
  {
    return (this->give_masterel() == fe) ? this->give_attributes()->give_cloads() : 0;
  }
  
};



//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS EDGE                          *** *** *** ***
//*  ********************************************************************************************
class Edge : public Facedge //, public ComponentGeometry1D
{
  // nno = nn               // number of nodes on Edge
  // ned = 0                // number of edges on Edge
  // nfa = 0                // number of faces on Edge
  // numsuperelem = NAEL_ED   // number of superior elements
  // numsuperface = NAFA_ED   // number of superior faces
  
  GPA<const Face>    superfaces;      // array of pointers to superior faces; this array should be sorted!!!
  
 public:
  //  Edge (long gid, const Mesh *mg, long nn) : Facedge (gid, gid, -1, mg, nn, 0, 0) { }
  Edge (long nn, const Point **nods) : Facedge (-1, -1, 0, NULL, nn, 0, 0) //, ComponentGeometry1D ()
  {
    superfaces.initialize(0,2);
    
    cg = new ComponentGeometry1D(this, ZERO);
    while (nn--)  points.assign(nn, nods[nn]);
    
    attributes = new FacedgeAttribs(this, (long)0);
  }
  //Edge () : ... () {    superelems.initialize(0,2); }
  virtual ~Edge() { }
  
  virtual classID give_classid() const { return classEdge; }
  
  //*  *** SET ***
  void setadd_superface (const Face    *comp, bool uniquecheck) { if (uniquecheck) superfaces.add_unique(comp); else superfaces.add(comp); } // ma to uniquecheck smyslo? to kdyz se dela re connectivity, ale na co to je??
  void remove_superface (const Face    *comp)                   {                  superfaces.free_ptr(comp);  superfaces.shake_down(); }
  
  //*  *** GET ***
  long                   give_numsuperface (void)   const { check_connectivity(); return superfaces(); }
  const Face*               give_superface (long i) const { check_connectivity(); return superfaces[i]; }
  const GPA<const Face   >* give_superfaces (void)  const { check_connectivity(); return &superfaces; }
  
  
  //*  *** CONNECTIVITY ***
  void connectivity_assembling (bool re=false);
  void connectivity_removing (void);
  void switch_myself_at_connectivity (Cell *master);
  
  virtual void set_model_prop (long val, const Model *model, bool flag=false);
  
 protected:
  virtual void integrate_duplicated_cell (const Edge *slave);
  
 public:
  void transform_nc (const Point* t1, const Point* t2, double &ksi);
  bool is_nod_on (double norm, const PoinT *point, double &ksi) const;
  const Face* give_superface (long ne, const Edge **edgs) const;
  virtual void switch_node_pointer (Point *slave, Point *master, bool duplcheck);
  
  virtual long give_ord (void) const { if (this->give_nno() != 2) errol; return -4; }
  
  double give_elemSize (void) const;
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const;
};

//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS FACE                          *** *** *** ***
//*  ********************************************************************************************
class Face : public Facedge //, public ComponentGeometry2D
{
  // nno = nn               // number of nodes on Face
  // ned = ne               // number of edges on Face
  // nfa = 0                // number of faces on Face
  // numsuperelem = 2         // number of superior elements
  // numsuperface = 0         // number of superior faces
  
  
 public:
  //Face (long gid, const Geometry *mg, long nn, long ne) : Facedge (gid, gid, -1, mg, 1, nn, ne, 0) { }
  Face (long nn, long ne, const Point **nods, const Edge **edgs) : Facedge (-1, -1, 0, NULL, nn, ne, 0) //, ComponentGeometry2D ()
  {
    if      (nn == 3)  cg = new ComponentGeometry2Dtriangle  (this, ZERO);
    else if (nn == 4)  cg = new ComponentGeometry2Dquadrangle(this, ZERO);
    else if (nn  > 4)  cg = new ComponentGeometry2Dpolygon   (this, ZERO);
    else  errol;
    
    while (nn--)  points.assign(nn, nods[nn]);
    while (ne--)  edges.assign(ne, edgs[ne]);
    
    attributes = new FacedgeAttribs(this, (long)0);
  }
  virtual ~Face() { }
  
  virtual classID give_classid() const { return classFace; }
  
  //*  *** CONNECTIVITY ***
  void connectivity_assembling (bool re=false);
  void connectivity_removing (void);
  void switch_myself_at_connectivity (Cell *master);
  
  virtual void set_model_prop (long val, const Model *model, bool flag=false);
  
  void transform_nc (const Edge *t1, const Edge *t2, PoinT *nc) const;
  virtual void switch_node_pointer (Point *slave, Point *master, bool duplcheck);
  
  virtual long give_ord (void) const { if (this->give_nno() != 3 && this->give_nno() != 4) errol; return -4; }
  
  double give_elemSize (void) const;
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const;
};


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS ELEMENT                         *** *** *** ***
//*  ********************************************************************************************
class Element : public Cell
{
 public:
  Element (classID mecg, long gid, long oid, long prop, const Geometry *mg, long nn, long ne, long nf) : Cell (mecg, gid, oid, prop, mg, nn, ne, nf) { }
  Element (classID mecg, const Element * src) : Cell (mecg, src) { }
  virtual ~Element() {}
  
 protected:
        ElemAttribs *elemAttribs()       { return (ElemAttribs*)attributes; }
  const ElemAttribs *elemAttribs() const { return (ElemAttribs*)attributes; }
  
 public:
  VTKPDtopology give_VTKPDtopology (void) const;
  const Facedge* give_same_dimension_facedge (void) const;
  //const Facedge* give_same_dimension_facedge_uniq (void) const;
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const = 0;
  virtual SStype      give_SSType_default      (void) const = 0;
  
  //*  *** SET ***
  void set_mprop (long val);
  void set_prop_node_inher (bool everynode);
  void set_fullhinge_at (const Point *p) { this->elemAttribs()->set_fullhinge_at(this->give_point_indx(p)); }
  
  //*  *** GET ***
        ElemAttribs *give_elemAttribs (void)       { return (ElemAttribs*)attributes; }
  const ElemAttribs *give_elemAttribs (void) const { return (ElemAttribs*)attributes; }
  bool has_fullhinge_at (int point) const;
  
  void connectivity_assembling (bool re=false);
  void connectivity_removing (void);
  
  virtual long give_edge_nodes       (const Point **&edgnodes) const { errol; return 0; }
  //
  void set_load  (int i, int indx=-1);
  
  
 protected:
  virtual void integrate_duplicated_cell (const Element *slave);
  
 public:
  bool is_point_on (const PoinT *point, const GeometryComponent *&comp, PoinT *nc) const;
  void switch_node_pointer_in_all_components (Point *slave, Point *master, bool duplcheck);
  virtual void switch_node_pointer (Point *slave, Point *master, bool duplcheck);
  
  //*  *** COMPUTE ***
  virtual double give_quality (void) const { return -1; }
  virtual double give_volume (void) const { errol; return 0.0; }
  
  //* tyhle dve fce se mi moc nelibi, nejsou systematicke, ale zatim nevim jak to udelat lepe
  bool isBeam  (void) const { return (this->give_dimension() == 1 && elemAttribs()->give_dpn() == DPN_DxyzRxyz) ? true : false; }
  bool isTruss (void) const { return (this->give_dimension() == 1 && elemAttribs()->give_dpn() == DPN_DxyzR___) ? true : false; }
  
  //*  *** PRINT ***
  void print_row_VTK (FILE *stream) const;
  virtual void print_row_VTX (char *str) const;
};



//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS GELEMENT                        *** *** *** ***
//*  ********************************************************************************************
class Gelement : public Element
{
 protected:
  GPA<const Vertex>   *fixedVertices;   // vertices fixed to gelement
  GPA<const Gelement> *fixedGelements;  // gelements fixed to gelement
  
 public:
  Gelement (classID mecg, long gid, long oid, const Geometry *mg, long nn, long ne, long nf);
  Gelement (classID mecg, const Gelement *src);
  virtual ~Gelement() { delete fixedVertices; delete fixedGelements; }
  
  //*  *** GET ***
        GelemAttribs *give_gelemAttribs (void)       { return (GelemAttribs*)attributes; }
  const GelemAttribs *give_gelemAttribs (void) const { return (GelemAttribs*)attributes; }
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const { return Pd->give_DOFspnod(); }
  virtual SStype      give_SSType_default      (void) const { return SST_Void; }
  
  const GPA<const Vertex>*   give_fixedVertices (void) const { return fixedVertices; }
  const GPA<const Gelement>* give_fixedGelements(void) const { return fixedGelements; }
  
  //
  void set_elemSize (double val);
  void set_elemCount (double val);
  
  //
  void assign_fixed_entities_by_ID (bool node, long ncn, const long *icn);
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const { errol; }
};


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS POLYLINE                        *** *** *** ***
//*  ********************************************************************************************
class PolyLine : public Gelement
{
 public:
  PolyLine (long gid, long oid, const Geometry *mg, int nv, classID mecg=classVoid) : Gelement (((mecg == classVoid) ? classComponentGeometry1Dpoly : mecg), gid, oid, mg, nv, (nv ? nv-1 : 0), 0) { }
  PolyLine (const PolyLine *src, classID mecg=classVoid) : Gelement (((mecg == classVoid) ? classComponentGeometry1Dpoly : mecg), src) { }
  virtual ~PolyLine() { }
  
  virtual classID give_classid() const { return classPolyLine; }
  virtual void checkConsistency (void) const;
  
  void add_point (const Point *val) {
    points.add(val);
    edges.resize(points()-1);
    ((ComponentGeometry1Dpoly*)cg)->reset_nv();
  }
  void del_points_down (long val)   {
    points.resizedown(val);
    edges.resizedown(points()-1);
    ((ComponentGeometry1Dpoly*)cg)->reset_nv();
  }
  
  virtual long give_edge_nodes (const Point **&edgnodes) const { errol; return 0; }
  
  void      init_point_on          (Mesh* pm, GPA<const Node> &hnodes, const PoinT *point, const GeometryComponent* comp, const PoinT *nc, const Point *parentpoint) const;
  PolyLine* generate_mesh_RFbyHN   (Mesh* pm, GPA<const Node> &hnodes) const;
  void      divide                 (Mesh* pm, GPA<const Node> &hnodes) const;
  void      findout_segment_domain (Mesh* pm, GPA<const Node> &hnodes) const;
  
  //*  *** PRINT ***
};

//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS LINE                          *** *** *** ***
//*  ********************************************************************************************
class Line : public PolyLine
{
 public:
  Line (long gid, long oid, const Geometry *mg) : PolyLine (gid, oid, mg, 2, classComponentGeometry1D) { }
  virtual ~Line() { }
  
  virtual classID give_classid() const { return classLine; }
  virtual void initialize (void);
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const { return 0; }
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const { }
};



//*  ********************************************************************************************
//*  *** *** *** ***                       CLASS POLYGONMDL                       *** *** *** ***
//*  ********************************************************************************************
class PolygonMdl : public Gelement //, public ComponentGeometry2Dtriangle
{
 public:
  PolygonMdl (long gid, long oid, const Geometry *mg) : Gelement (classComponentGeometry2Dpolygon, gid, oid, mg, 1, 1, 1) { }
  virtual ~PolygonMdl() { }
  
  virtual classID give_classid() const { return classPolygonMdl; }
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const;
  
  void allocNnNe (int n);
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const { }
};


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS FELEMENT                        *** *** *** ***
//*  ********************************************************************************************
class FElement : public Element
{
 protected:
  long domain;                   // domain
  long lid;                      // local identification number of element
  
  const Gelement *mdl_masterel;  // master element from master model
  
  
 public:
  FElement (classID mecg, long gid, long oid, const Geometry *mg, long ord, long nn, long ne, long nf, bool aa, long dom, long lid);
  FElement (classID mecg, const FElement * src);
  virtual ~FElement();
  
  void attributes_allocation (const GelemAttribs *masterat);
  
  virtual void initialize (void);
  virtual void checkConsistency (void) const;
  
  //
  const Mesh* Msh(void) const { return (Mesh*)Geom; }
  
  virtual void set_model_prop (long val, const Model *model, bool flag=false);
  
  //virtual FiniteElementType give_fe (enumSolver sol) const = 0;
  virtual SPRpatchType give_spr_patch_type (void) const { errol; return SPRPT_Void; }
  int give_nDOFs (void) const { return this->give_nno() * DOFsPerNode2nDOFs(this->elemAttribs()->give_dpn()); }
  IntPointSet give_IPset_comp (Solver sol) const;
  IntPointSet give_IPset_rslts (Solver sol) const;
  //int         give_NIP   (enumSolver sol, ResultTypesAtElem rte) const { return IntPointSet_e2i (this->give_IPset(sol), sol, rte); }
  int         give_rslt_NIP (Solver sol) const { return IntPointSet_give_number_ips (this->give_IPset_rslts(sol)); }
  //IntPointSet_give_number_ips_result (, sol, RTE_stress); }
  
  IntPointSet give_IPset2 (void) const;
  
  //*  *** SET ***
  void set_lid    (long val) { lid            = val; }
  int  give_ord   (void) const { return this->elemAttribs()->give_ord(); }
  long give_nnoed (int i) const { if (this->elemAttribs()->give_ord() != -4)  errol;  return 2; } // number of nodes on edge
  //*  *** GET ***
  long give_domain (void) const { return domain; }
  long give_lid    (void) const { return lid; }
  //bool is_pl_orig  (void) const { return ( (this->mdl_masterel && this->mdl_masterel->give_Geom()->give_gp() == MMP_RFbyHN)); }
  const Gelement* give_mdl_masterel (void) const { return mdl_masterel; }
  
  //*  *** SET ***
  void set_node (long i,       long  nid) { points.assign(i, (Point*)Msh()->give_Node(nid)); }
  void set_node (long i, const Node* nod) { points.assign(i, (Point*)nod); }
  //*  *** GET ***
  const Node* give_node (long i) const {                       return (const Node*)points[i]; }
  //const GPA<const Node>* give_nodes (void) const {                       return &nodes; } // !!!! zkusit vrati &
  // pokud uz existuje "master geometry", tak existuji i mdl_masterel-y elementu
  //long give_parent_id (void) const { _errorr("tuto fci opravit"); return (mdl_masterel ?  : this->give_origID() ); }
  long give_parent_id   (void) const { if (!mdl_masterel) _errorr("existuje master model?");  return this->give_mproperty() - 1; } // mdl_masterel->give_origID()
  long give_parent_prop (void) const { if (!mdl_masterel) _errorr("existuje master model?");  return mdl_masterel->give_mpropertyORzero(); }
  
 protected:
  virtual void integrate_duplicated_cell (const Element *slave);
  
  
  //*  *******************
  //*  ***   RESULTS   ***
  //*  *******************
 private:
  Array ***results;
  long regid;          // adaptivity region number
 protected:
  double *maxSigmaEq;  // maximal stress equivalent at element
  
 public:
  //*  *** SET ***
  void set_regid (long val) { regid = val; }
  //void set_resultN (long s, double *rslt,             long step, ResultTypesAtNode rt);


  //*  *** GET ***
  long give_regid (void) const { return regid; }
  void allocate_results (void);
  void add_result (Array *rslt, long step, ResultTypesAtElem rt);
  void set_result (long s, double *rslt, long step, ResultTypesAtElem rt);
  void set_result (        double  rslt, long step, ResultTypesAtElem rt);
  void set_result (  const VectoR *rslt, long step, ResultTypesAtElem rt);
  const Dscal* give_results_ds (long step, ResultTypesAtElem rt) const;
  const Dvctr* give_results_dv (long step, ResultTypesAtElem rt) const;
  const Dmtrx* give_results_dm (long step, ResultTypesAtElem rt) const;
  const Array* give_results    (long step, ResultTypesAtElem rt) const { return results[step][rt]; }
  // //double give_displ (int i, int j) const { return displ[i][j]; }
  // Dvctr* give_full_alloc_displcmnt_at (long step) const;
  // Dvctr* give_full_alloc_reactions_at (long step) const;
  // void give_displcmnt_in_all_steps(Dvctr *data, int indx) const;
  virtual void compute_maxSigmaEq (void) { _errorr2("CS usage computation is not supported at element %ld", ID+1); }
  void setup_maxSigmaEq (void);
  double give_CSusage_elast (void);
  double give_CSusage_elast_rel (void);
  int    give_CSusage_elast_bool (void);
  
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL) = 0;
  //const Dmtrx** give_plaststrain (void) const { return (const Dmtrx**)plaststrain; }
  //virtual void rotate2local(Dvctr *data) { _errorr ("not implemented"); }
  void check_rslts (long step, ResultTypesAtElem rt) const { if (!results || !results[step][rt])
      _errorr3 ("no result %d at step %ld", (int)rt, step); }
  double fillupbyzero (Dvctr *data, SStype SST) const;
  long give_result_ncomp (long time_step, ResultTypesAtElem rte) const;
  
  
  void initialize_CSL (int c, long *fID, long *tCSL);
  
  //*  *** READ ***
  virtual void read_input (const char *&str, femFileFormat fff);
  virtual void read_output_OOFEM (FILE *stream, long step)                       { _errorr ("not implemented"); }
  virtual void read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt) { _errorr ("not implemented"); }
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const;
  //void print_row_VTK (FILE *stream) const;
  //virtual void print_row_VTX (char *str) const;
  
  //*  *** READ ***
  void read_nodes (const char *&str, femFileFormat fff);
  //*  *** PRINT ***
  
  bool is_cmfr (void);
  
  virtual void give_ip_coords_global (IntPointSet ips, int i, PoinT &coords) const { _errorr("not implemented"); }
};


//*  ********************************************************************************************
//*  *** *** *** ***                          CLASS BEAM                          *** *** *** ***
//*  ********************************************************************************************
class Beam : public FElement //, public ComponentGeometry1D
{
 public:
  Beam (long gid, long oid, const Geometry *mg, bool aa, long dom, long lid) : FElement (classComponentGeometry1D, gid,oid, mg,1  ,2 ,1 ,0 ,aa,dom,lid) {}
  virtual ~Beam() { }
  
  virtual classID give_classid() const { return classBeam; }
  virtual void checkConsistency (void) const;
  
  virtual void initialize (void) {
    FElement :: initialize ();
  }
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const;
  virtual SStype      give_SSType_default      (void) const;
  //virtual FiniteElementType give_fe (Solver sol) const;
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const { return 0; }
  void give_vector (VectoR *v) const;
  
  
  //*  *******************
  //*  ***   RESULTS   ***
  //*  *******************
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL);
  double give_normal_force (long step) const;
  virtual void compute_maxSigmaEq (void);
  
  //*  *** READ ***
  void read_output_OOFEM (FILE *stream, long step);
  void read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt);
  
  //*  *** PRINT ***
  virtual void print_row (FILE *stream, femFileFormat fff, bool endline=true, long did=0) const;
};


//*  ********************************************************************************************
//*  *** *** *** ***                        CLASS TRIANGLE                        *** *** *** ***
//*  ********************************************************************************************
class Triangle : public FElement //, public ComponentGeometry2Dtriangle
{
 protected:
  Dmtrx *rotMg2l;  // rotation matrix from global to local SS
  
 public:
  Triangle (long gid, long oid, const Geometry *mg, bool aa, long dom, long lid) : FElement (classComponentGeometry2Dtriangle, gid,oid,mg,1  ,3 ,3 ,1 , aa, dom,lid) { rotMg2l = NULL; }
  Triangle (const Quadrangle * src, bool firstr, bool firstdiag);
  virtual ~Triangle() { delete rotMg2l; }
  
  virtual classID give_classid() const { return classTriangle; }
  virtual void checkConsistency (void) const;
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const;
  virtual SStype      give_SSType_default      (void) const;
  //virtual FiniteElementType give_fe (Solver sol) const;
  virtual SPRpatchType give_spr_patch_type (void) const { return SPRPT_2d; }
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const;
  
  //*  *******************
  //*  ***   RESULTS   ***
  //*  *******************
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL);
  virtual void compute_maxSigmaEq (void);
  
  //*  *** COMPUTE ***
  virtual double give_quality (void) const;
  virtual double give_volume (void) const { return this->cg->give_lav() * this->elemAttribs()->give_cs()->give_thickness(); }
  
  void compute_rotMg2l(void);
  virtual void rotate2local(Dvctr *data);
  
  virtual void give_ip_coords_global (IntPointSet ips, int i, PoinT &coords) const;
  
  //*  *** READ ***
  void read_input (const char *&str, femFileFormat fff);
  void read_output_OOFEM (FILE *stream, long step);
  void read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt);
};


//*  ********************************************************************************************
//*  *** *** *** ***                      CLASS QUADRANGLE                      *** *** *** ***
//*  ********************************************************************************************
class Quadrangle : public FElement //, public ComponentGeometry2Dquadrangle
{
 protected:
  Dmtrx *rotMg2l;  // rotation matrix from global to local SS
  
 public:
  Quadrangle (long gid, long oid, const Geometry *mg, bool aa, long dom, long lid) : FElement (classComponentGeometry2Dquadrangle, gid,oid,mg,1  ,4 ,4 ,1 , aa, dom,lid)
  {
    rotMg2l = NULL;
  }
  virtual ~Quadrangle() { delete rotMg2l; }
  
  virtual classID give_classid() const { return classQuadrangle; }
  virtual void checkConsistency (void) const;
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const;
  virtual SStype      give_SSType_default      (void) const;
  //virtual FiniteElementType give_fe (Solver sol) const;
  virtual SPRpatchType give_spr_patch_type (void) const { return SPRPT_2d; }
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const;
  
  bool is_first_diag_short (void);
  
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL);
  
  //*  *** COMPUTE ***
  //virtual double give_quality (void);
  virtual double give_volume (void) const { return this->cg->give_lav() * this->elemAttribs()->give_cs()->give_thickness(); }
  
  //void compute_rotMg2l(void);
  //virtual void rotate2local(Dvctr *data);
  
  virtual void give_ip_coords_global (IntPointSet ips, int i, PoinT &coords) const;
  
  void read_input (const char *&str, femFileFormat fff);
  void read_output_OOFEM (FILE *stream, long step);
};


//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS TETRA                          *** *** *** ***
//*  ********************************************************************************************
class Tetra : public FElement
{
  long *surfprop;  // property of surfaces (specially for t3d mesh input)
  
 public:
  Tetra (long gid, long oid, const Geometry *mg, bool aa, long dom, long lid) : FElement (classComponentGeometry3Dtetrahedron, gid,oid,mg,1  ,4 ,6 ,4 , aa, dom,lid)
  {
    surfprop = NULL;
  }
  virtual ~Tetra() { delete [] surfprop; }
  
  virtual classID give_classid() const { return classTetra; }
  virtual void checkConsistency (void) const;
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const;
  virtual SStype      give_SSType_default      (void) const;
  //virtual SStype give_SStype (void) const {
  //  if (attributes->give_sst() == SST_3d)  return RSLTSST_3dshell_nodes;
  //  else errol;  return RSLTSST_Void;
  //}
  //virtual FiniteElementType give_fe (Solver sol) const;
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const;
  
  //void renumbering(void);
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL) { errol; return 0.; }
  
  //*  *** READ ***
  //void read_output_OOFEM (FILE *stream, long step);
  //virtual void read_T3d_input (const char *str);
  
  //*  *** PRINT ***
  virtual void print_row_VTX (char *str) const;
};

//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS BRICK                          *** *** *** ***
//*  ********************************************************************************************
class Brick : public FElement //, public ComponentGeometry3Dhexahedron
{
 public:
  Brick (long gid, long oid, const Geometry *mg, bool aa, long dom, long lid) : FElement (classComponentGeometry3Dhexahedron, gid,oid,mg,1  ,8 ,12,6 , aa, dom,lid) //, ComponentGeometry3Dhexahedron ()
  { }
  virtual ~Brick() { }
  
  virtual classID give_classid() const { return classBrick; }
  virtual void checkConsistency (void) const;
  
  virtual DOFsPerNode give_DOFsPerNode_default (void) const;
  virtual SStype      give_SSType_default      (void) const;
  //virtual SStype give_SStype (void) const {
  //  if (attributes->give_sst() == SST_3d)  return RSLTSST_3dshell_nodes;
  //  else errol;  return RSLTSST_Void;
  //}
  //virtual FiniteElementType give_fe (Solver sol) const;
  
  virtual long give_edge_nodes (const Point **&edgnodes) const;
  virtual long give_face_nodes_edges (const Point **&facnodes, const Edge **&facedges) const;
  
  void renumbering(void);
  virtual double give_ssstate (Dvctr *data, SStype SST, RVType rvtype, char type, long step, const Node *node = NULL);
  
  //*  *** READ ***
  void read_output_OOFEM (FILE *stream, long step);
  void read_output_SIFEL (FILE *stream, long step, ResultTypesAtElem rt);
  virtual void read_input (const char *&str, femFileFormat fff);
  
  //*  *** PRINT ***
  virtual void print_row_VTX (char *str) const;
};


} // namespace midaspace

#endif // MIDAS_CELL_H