src/tetmesh/tet.cpp

Go to the documentation of this file.

/***************************************************************************
                          tet.cpp  -  description
                             -------------------
    begin                : Fri Dec 12 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <stdexcept>
#include "vector.h"
#include "linalg.h"
#include "point.h"
#include "tet.h"
#include "tetmesh.h"
#include "fmxxtopology.h"

namespace mesh {

using namespace colarray;

  Tet::Tet(const TetMesh& mesh, id_t id, id_t corners[4], id_t material) :
    Entity(id), _edge_orientation(0), _material(material)
{
    // store corners
    for (int i = 0; i < 4; i ++)
        _corners[i] = corners[i];

    // calculate Tet volume
    double volume = get_volume();
    
    // check for degenerate elements
    if (1.0 + volume == 1.0)
        throw std::runtime_error("Tetrahedron has zero volume.");

    // enforce positive orientation
    if (volume < 0.0)
        std::swap<id_t>(_corners[0], _corners[1]);
}

void Tet::set_face_id(int i, id_t id) {
    _faces[i] = id;
}

void Tet::set_edge_id(int i, id_t id) {
    _edges[i] = id;
}

void Tet::set_edge_orientation(int edge_lid, bool orientation) {
    if (orientation)
        _edge_orientation |= 1 << edge_lid;
    else
        _edge_orientation &= ~(1 << edge_lid);
}

id_t Tet::get_material_ID() const { 
  return _material;
}

double Tet::get_volume() const { 
    CoordinateMatrix X(*this);
    return linalg::determinant(X) / 6.0;
}

Vector3 Tet::get_center_of_gravity() const
{
        int i; 
        Vector3 cog(0,0,0);
        
        for(i=0;i<NCORNERTET;i++)
        {
                cog+=0.25 * this->get_corner(i)->get_coord();
        }
                
        /* go back */
        return(cog);
}


void Tet::cartesian_to_simplex(const Vector3& cartesian, Vector4& simplex) const {
    CoordinateMatrix X(*this);
    ColumnVector<double> t(3);

    for (int i = 0; i < 4; i ++) {
        RowVector<double> r(X, i);           // create alias to row i
        Vector<double> r1(r, 1, 4);          // skip 1st element in row i
        t.inject(r1);                        // save row i
        r1(0) = cartesian.x;                 // overwrite row i with x
        r1(1) = cartesian.y;
        r1(2) = cartesian.z;
        simplex[i] = linalg::determinant(X) / (6.0 * get_volume()); // _volume;
        r1.inject(t);                        // restore row i
    }
}

bool Tet::is_inside(const Vector3& point) const {
    Vector4 simplex_coord;
    cartesian_to_simplex(point, simplex_coord);
    for (int i = 0; i < 4; i ++)
        if (simplex_coord[i] < 0.0)
            return false;
    for (int i = 0; i < 4; i ++)
        if (simplex_coord[i] > 1.0)
            return false;
    return true;
}

bool Tet::is_inside_gracious(const Vector3& point) const {
    const double grace_tolerance = 1e-14;
    Vector4 simplex_coord;

    cartesian_to_simplex(point, simplex_coord);
    for (int i = 0; i < 4; i ++)
        if (simplex_coord[i] < -grace_tolerance)
            return false;
    for (int i = 0; i < 4; i ++)
        if (simplex_coord[i] > 1.0 + grace_tolerance)
            return false;
    return true;    
}

void Tet::get_bounding_box(Vector3& center, Vector3& extent) const {
    Vector3 min, max, coord;

    min = get_corner(0)->get_coord();
    max = min;
    for (int i = 1; i < 4; i ++) {
        coord = get_corner(i)->get_coord();
        min.make_floor(coord);
        max.make_ceil(coord);
    }
    center = max.mid_point(min);
    extent = 0.5 * (max - min);
}

Tet::CoordinateMatrix::CoordinateMatrix(const Tet& tet) : colarray::Matrix<double>(4, 4) {
    Vector3 coord;

    for (int i = 0; i < 4; i ++) {
        this->operator()(i,0) = 1.0;
        coord = tet.get_mesh()->get_point(tet._corners[i])->get_coord();
        for (int j = 0; j < 3; j ++)
            this->operator()(i,j+1) = coord[j];
    }
}

Tet::GradientMatrix::GradientMatrix(const CoordinateMatrix& X) : colarray::Matrix<double>(3, 4) {
    colarray::Matrix<double> Xinv(4, 4);
    // compute inverse of coordinate matrix
    linalg::inverse(X, Xinv);
    // rows 1, 2, 3 are gradients
    this->inject(colarray::Matrix<double>(Xinv, 1, 4, 0, 4));
}

Tet::GradientMatrix::GradientMatrix(const Tet& tet) : colarray::Matrix<double>(3, 4) {
    Tet::CoordinateMatrix X(tet);
    colarray::Matrix<double> Xinv(4, 4);
    // compute inverse of coordinate matrix
    linalg::inverse(X, Xinv);
    // rows 1, 2, 3 are gradients
    this->inject(colarray::Matrix<double>(Xinv, 1, 4, 0, 4));
}

Point* Tet::get_corner(int i) const {
    return get_mesh()->get_point(get_corner_id(i));
}

Edge* Tet::get_edge(int i) const {
    return get_mesh()->get_edge(get_edge_id(i));
}

Face* Tet::get_face(int i) const {
    return get_mesh()->get_face(get_face_id(i));
}

TetMesh* Tet::get_mesh() const {
    return TetMesh::get_instance();
}

} // namespace mesh

---