compgeom.cpp

Go to the documentation of this file.

#include "compgeom.h"

#include "point.h"

#include "geometrylib.h"

#include <math.h>


namespace midaspace {

void ComponentGeometry :: initialize (void)
{
  initialized = true;
  
  const GPA<const Point> *points = owner->give_points();
  
  int i;
  
  // centrcoord
  centrcoord.copy(points[0][0]->give_coords());
  for (i=1; i<points[0](); i++)  centrcoord.add(points[0][i]->give_coords());
  centrcoord.dvd(points[0]());
  
  double d;
  // circum
  circum = 0.0;
  for (i=0; i<points[0](); i++) {
    d = centrcoord.dist2_to(points[0][i]->give_coords());
    if (circum < d)
      circum = d;
  }
  circum = (1.0 + ZERO) * sqrt(circum);
}
  

void ComponentGeometry1D :: initialize (void)
{
  ComponentGeometry :: initialize();
  
  const GPA<const Point> *points = owner->give_points();
  
  if (points[0]() == 2)  lav = points[0][0]->give_coords()->dist_to(points[0][1]->give_coords());
  else                   _errorr ("Unsupported number of nodes on 1d element.");
}

void ComponentGeometry1Dpoly :: initialize (void)
{
  ComponentGeometry :: initialize();
  
  const GPA<const Point> *points = owner->give_points();
  
  lav = 0.0;
  
  for (int i=1; i<nv; i++)
    lav += points[0][i-1]->give_coords()->dist_to(points[0][i]->give_coords());
}

ComponentGeometry1Dpoly :: ComponentGeometry1Dpoly (Cell *o, double Z) : ComponentGeometry1D(o, Z)
{ nv = owner->give_nno(); }

void ComponentGeometry1Dpoly :: checkConsistency (void) const
{
  ComponentGeometry1D :: checkConsistency();
  if (owner->give_nno() != nv)  errol;
}

void ComponentGeometry1Dpoly :: reset_nv (void) { nv = owner->give_nno(); }


ComponentGeometry2D :: ComponentGeometry2D (Cell *o, double Z) : ComponentGeometry(o, Z)
{ nv = owner->give_nno(); }

void ComponentGeometry2D :: checkConsistency (void) const
{
  ComponentGeometry :: checkConsistency();
  if (owner->give_nno() != nv)  errol;
}

void ComponentGeometry2D :: initialize (void)
{
  ComponentGeometry :: initialize();
  
  const GPA<const Point> *points = owner->give_points();
  
  // normal vector
  if (nv == 4)  normal.beVectProduct (points[0][0]->give_coords(), points[0][2]->give_coords(), points[0][1]->give_coords(), points[0][3]->give_coords());
  else          normal.beVectProduct (points[0][0]->give_coords(), points[0][1]->give_coords(), points[0][2]->give_coords());
  
  // check collapse
  double l = normal.give_length();
  if (!owner->give_Pd()->give_PDBO(PDBO_MultipElem)  &&  l < circum * circum * zero)
    _warningg2("collapsed 2d element %ld", owner->give_ID()+1);
  
  double d = -1 * normal.giveScalProduct (points[0][0]->give_coords());
  double dist;
  
  //double zero_prec = zero;     // 
  //double zero_prec = 1e-8;       // docasne, kvuli presnosti t3d
  for (int i=3; i<nv; i++) {
    dist = (normal.giveScalProduct(points[0][i]->give_coords()) + d) / l;
    if (dist < -1 * circum * zero  ||  circum * zero < dist)
      _errorr3("point %ld is not in plane of element %ld", i+1, owner->give_ID()+1);
  }
  
  if (nv <= 4)  lav = 0.5 * l;
  else {
    const PoinT **pointss = new const PoinT*[nv+2];
    
    for (int i=0; i<nv; i++)  pointss[i] = points[0][i]->give_coords();
    pointss[nv  ] = points[0][0]->give_coords();
    pointss[nv+1] = points[0][1]->give_coords();
    
    VectoR norm;
    norm.copy(&normal)->normalize();
    
    lav =  area3D_Polygon (nv, pointss, &norm);
    
    // oprava smeru normaly u konkavniho polygonu
    const Element *elem = dynamic_cast<Element*>(owner);
    if (!elem) {
      Face *face = dynamic_cast<Face*>(owner);
      if (!face) errol;
      elem = face->give_main_masterel_uniq();
    }
    if (elem) {
      const LCS_beam *llccss = elem->give_elemAttribs()->give_lcs_global();
      if (llccss) {
    const Elem3D* vect = llccss->give_vector();
    if (!norm.is_identical_to (vect, 0.000001)) { norm.tms(-1.0); normal.tms(-1.0); }
    if (!norm.is_identical_to (vect, 0.000001)) 
      errol;
      }
    }
    else errol;
    
    delete [] pointss;
  }
}

const LCS_pure* ComponentGeometry2D :: give_lcs (void)
{
  if (lcs)  return lcs;
  
  //
  const GPA<const Point> *points = owner->give_points();
  
  // normal vector
  //lcs->v1.beVectProduct (points[0][0]->give_coords(), points[0][1]->give_coords(), points[0][2]->give_coords());
  VectoR v1, v2, v3;
  
  v1.beP2P( points[0][0]->give_coords(), points[0][1]->give_coords()  );
  v2.beP2P( points[0][1]->give_coords(), points[0][2]->give_coords()  );
  
  v3.beVectProduct( &v1, &v2 );
  v2.beVectProduct( &v3, &v1 );
  
  v1.normalize();
  v2.normalize();
  v3.normalize();
  
  lcs = new LCS_pure;
  lcs->set_n1(&v1);
  lcs->set_n2(&v2);
  lcs->set_n3(&v3);
  
  return lcs;
}



void ComponentGeometry3Dtetrahedron :: initialize (void)
{
  ComponentGeometry :: initialize();
  
  const GPA<const Point> *points = owner->give_points();
  
  const PoinT *c1 = points[0][0]->give_coords();
  const PoinT *c2 = points[0][1]->give_coords();
  const PoinT *c3 = points[0][2]->give_coords();
  const PoinT *c4 = points[0][3]->give_coords();
  
  Dmtrx m(3,3);
  
  m(0,0) = c1->x - c4->x;   m(0,1) = c2->x - c4->x;   m(0,2) = c3->x - c4->x;
  m(1,0) = c1->y - c4->y;   m(1,1) = c2->y - c4->y;   m(1,2) = c3->y - c4->y;
  m(2,0) = c1->z - c4->z;   m(2,1) = c2->z - c4->z;   m(2,2) = c3->z - c4->z;
  
  lav = m.give_determinant() / 6.0;
  if (lav < 0.0)  errol;
}
void ComponentGeometry3Dhexahedron :: initialize (void)
{
  ComponentGeometry :: initialize();
  
  const GPA<const Point> *points = owner->give_points();
  
  int i;
  Dvctr x(8),y(8),z(8);
  
  // 
  for (i=0; i<8; i++) {
    x[i] = points[0][i]->give_coord(0);
    y[i] = points[0][i]->give_coord(1);
    z[i] = points[0][i]->give_coord(2);
  }
  
  // jac je zatim na sifel cislovani uzlu !!!!!!!!!
  jac_3d (lav, x, y, z, 0,0,0);
  
  lav *= -8.0;
}

} // namespace midaspace