582 lines
18 KiB
C++
582 lines
18 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _GCTL_TETRAHEDRON_H
|
|
#define _GCTL_TETRAHEDRON_H
|
|
|
|
#include "vertex.h"
|
|
#include "entity.h"
|
|
|
|
namespace gctl
|
|
{
|
|
// Declaration of the basic tetrahedron type
|
|
template <typename A> struct type_tetrahedron;
|
|
|
|
typedef type_tetrahedron<void> tetrahedron; // tetrahedron type of attribute type of void
|
|
|
|
static int regular_order[12] = {0,1,2,0,2,3,0,3,1,1,3,2};
|
|
static int inverse_order[12] = {0,1,3,0,2,1,0,3,2,1,2,3};
|
|
|
|
/**
|
|
* Structure of a tetrahedron
|
|
*
|
|
* regular type of tetrahedron
|
|
* 3
|
|
* /|\
|
|
* / | \ y
|
|
* z / | \ /
|
|
* | / 1 \
|
|
* | / / \ \
|
|
* |/ / \ \
|
|
* 0-------------2----> x
|
|
*
|
|
* triangle list (anti-clockwise)
|
|
* 0 1 2
|
|
* 0 2 3
|
|
* 0 3 1
|
|
* 1 3 2
|
|
*
|
|
* inverse type of tetrahedron
|
|
* 3
|
|
* /|\
|
|
* / | \ y
|
|
* z / | \ /
|
|
* | / 2 \
|
|
* | / / \ \
|
|
* |/ / \ \
|
|
* 0-------------1----> x
|
|
*
|
|
* triangle list (anti-clockwise)
|
|
* 0 1 3
|
|
* 0 2 1
|
|
* 0 3 2
|
|
* 1 2 3
|
|
*
|
|
* static int regular_order[12] = {0,1,2,0,2,3,0,3,1,1,3,2};
|
|
* static int inverse_order[12] = {0,1,3,0,2,1,0,3,2,1,2,3};
|
|
*/
|
|
template <typename A>
|
|
struct type_tetrahedron : public entity<vertex3dc, 4, A>
|
|
{
|
|
int *vec_order; ///< index of the local nodes anti-clockwise
|
|
type_tetrahedron<A> *neigh[4]; ///< index of the neighboring tetrahedrons
|
|
/**
|
|
* constructor
|
|
*/
|
|
type_tetrahedron();
|
|
/**
|
|
* @brief Set object from parameters
|
|
*
|
|
* @warning This function will locate memories to store vertice
|
|
*
|
|
* @param[in] p0 The first point
|
|
* @param[in] p1 The second point
|
|
* @param[in] p2 The third point
|
|
* @param[in] p3 The fourth point
|
|
* @param[in] index The element index
|
|
*/
|
|
type_tetrahedron(const point3dc &p0, const point3dc &p1, const point3dc &p2,
|
|
const point3dc &p3, int index = 0);
|
|
/**
|
|
* @brief Set object from parameters
|
|
*
|
|
* @warning This function will locate memories to store vertice
|
|
*
|
|
* @param[in] ps0 The first point
|
|
* @param[in] ps1 The second point
|
|
* @param[in] ps2 The third point
|
|
* @param[in] ps3 The fourth point
|
|
* @param[in] index The element index
|
|
*/
|
|
type_tetrahedron(const point3ds &ps0, const point3ds &ps1, const point3ds &ps2,
|
|
const point3ds &ps3, int index = 0);
|
|
/**
|
|
* @brief Constructor of the tetrahedron with initial parameters
|
|
*
|
|
* @param[in] index Index of the element
|
|
* @param vert0 The vertex 0
|
|
* @param vert1 The vertex 1
|
|
* @param vert2 The vertex 2
|
|
* @param vert3 The vertex 3
|
|
*/
|
|
type_tetrahedron(vertex3dc &vert0, vertex3dc &vert1, vertex3dc &vert2,
|
|
vertex3dc &vert3, int index = 0);
|
|
/**
|
|
* @brief Constructor of the tetrahedron with initial parameters
|
|
*
|
|
* @param[in] index Index of the element
|
|
* @param vert0 The vertex 0
|
|
* @param vert1 The vertex 1
|
|
* @param vert2 The vertex 2
|
|
* @param vert3 The vertex 3
|
|
*/
|
|
type_tetrahedron(vertex3dc *vertp0, vertex3dc *vertp1, vertex3dc *vertp2,
|
|
vertex3dc *vertp3, int index = 0);
|
|
/**
|
|
* @brief de-constructor
|
|
*/
|
|
virtual ~type_tetrahedron(){}
|
|
/**
|
|
* @brief Set object from parameters
|
|
*
|
|
* @warning This function will locate memories to store vertice
|
|
*
|
|
* @param[in] p0 The first point
|
|
* @param[in] p1 The second point
|
|
* @param[in] p2 The third point
|
|
* @param[in] p3 The fourth point
|
|
* @param[in] index The element index
|
|
*/
|
|
void set(const point3dc &p0, const point3dc &p1, const point3dc &p2,
|
|
const point3dc &p3, int index = 0);
|
|
/**
|
|
* @brief Set object from parameters
|
|
*
|
|
* @warning This function will locate memories to store vertice
|
|
*
|
|
* @param[in] ps0 The first point
|
|
* @param[in] ps1 The second point
|
|
* @param[in] ps2 The third point
|
|
* @param[in] ps3 The fourth point
|
|
* @param[in] index The element index
|
|
*/
|
|
void set(const point3ds &ps0, const point3ds &ps1, const point3ds &ps2,
|
|
const point3ds &ps3, int index = 0);
|
|
/**
|
|
* @brief Set the tetrahedron with input parameters
|
|
*
|
|
* @param[in] index Index of the element
|
|
* @param vert0 The vertex 0
|
|
* @param vert1 The vertex 1
|
|
* @param vert2 The vertex 2
|
|
* @param vert3 The vertex 3
|
|
*/
|
|
void set(vertex3dc &vert0, vertex3dc &vert1, vertex3dc &vert2,
|
|
vertex3dc &vert3, int index = 0);
|
|
/**
|
|
* @brief Set the tetrahedron with input parameters
|
|
*
|
|
* @param[in] index Index of the element
|
|
* @param vert0 The vertex 0
|
|
* @param vert1 The vertex 1
|
|
* @param vert2 The vertex 2
|
|
* @param vert3 The vertex 3
|
|
*/
|
|
void set(vertex3dc *vertp0, vertex3dc *vertp1, vertex3dc *vertp2,
|
|
vertex3dc *vertp3, int index = 0);
|
|
/**
|
|
* @brief Reset the structure to initial state
|
|
*/
|
|
void reset();
|
|
/**
|
|
* @brief Set neighbors of the tetrahedron
|
|
*
|
|
* @param nei_ptr0 The neighbor 0
|
|
* @param nei_ptr1 The neighbor 1
|
|
* @param nei_ptr2 The neighbor 2
|
|
* @param nei_ptr3 The neighbor 3
|
|
*/
|
|
void set_neighbor(type_tetrahedron<A> &nei0, type_tetrahedron<A> &nei1,
|
|
type_tetrahedron<A> &nei2, type_tetrahedron<A> &nei3);
|
|
/**
|
|
* @brief Get pointer of the j-th vertex on the i-th facet
|
|
*
|
|
* @param[in] i facet index (smaller than 4)
|
|
* @param[in] j vertex index (smaller than 3)
|
|
*
|
|
* @return vertex pointer
|
|
*/
|
|
vertex3dc *get(unsigned int i, unsigned int j) const;
|
|
/**
|
|
* @brief Get pointer of the j-th vertex on the i-th facet. Without any checks
|
|
*
|
|
* @param[in] i facet index (smaller than 4)
|
|
* @param[in] j vertex index (smaller than 3)
|
|
*
|
|
* @return { description_of_the_return_value }
|
|
*/
|
|
vertex3dc *fget(unsigned int i, unsigned int j) const;
|
|
/**
|
|
* @brief Initializes vec_order
|
|
*/
|
|
void deter_vert_order();
|
|
/**
|
|
* @brief Determines whether the specified tetrahedron is adjoined.
|
|
*
|
|
* @param nei_ptr The tetrahedron
|
|
*
|
|
* @return True if the specified tetrahedron is adjoined, False otherwise.
|
|
*/
|
|
bool is_adjoined(type_tetrahedron<A> &nei);
|
|
/**
|
|
* @brief Calculate the tetrahedron's volume
|
|
*
|
|
* @return volume
|
|
*/
|
|
double volume();
|
|
/**
|
|
* @brief Calculate the center of the tetrahedron
|
|
*
|
|
* @return point3dc center position
|
|
*/
|
|
point3dc center();
|
|
};
|
|
|
|
template <typename A>
|
|
type_tetrahedron<A>::type_tetrahedron() : entity<vertex3dc, 4, A>::entity()
|
|
{
|
|
vec_order = nullptr;
|
|
neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
|
|
}
|
|
|
|
template <typename A>
|
|
type_tetrahedron<A>::type_tetrahedron(const point3dc &p0, const point3dc &p1,
|
|
const point3dc &p2, const point3dc &p3, int index) : type_tetrahedron()
|
|
{
|
|
set(p0, p1, p2, p3, index);
|
|
}
|
|
|
|
template <typename A>
|
|
type_tetrahedron<A>::type_tetrahedron(const point3ds &ps0, const point3ds &ps1,
|
|
const point3ds &ps2, const point3ds &ps3, int index) : type_tetrahedron()
|
|
{
|
|
set(ps0, ps1, ps2, ps3, index);
|
|
}
|
|
|
|
template <typename A>
|
|
type_tetrahedron<A>::type_tetrahedron(vertex3dc &vert0, vertex3dc &vert1,
|
|
vertex3dc &vert2, vertex3dc &vert3, int index) : type_tetrahedron()
|
|
{
|
|
set(vert0, vert1, vert2, vert3, index);
|
|
}
|
|
|
|
template <typename A>
|
|
type_tetrahedron<A>::type_tetrahedron(vertex3dc *vertp0, vertex3dc *vertp1,
|
|
vertex3dc *vertp2, vertex3dc *vertp3, int index) : type_tetrahedron()
|
|
{
|
|
set(vertp0, vertp1, vertp2, vertp3, index);
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::set(const point3dc &p0, const point3dc &p1,
|
|
const point3dc &p2, const point3dc &p3, int index)
|
|
{
|
|
if (index < 0)
|
|
{
|
|
throw out_of_range("Invalid index number, From type_tetrahedron::set(...)");
|
|
}
|
|
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
this->vert[i] = new vertex3dc;
|
|
}
|
|
this->self_host = true;
|
|
|
|
this->id = index;
|
|
this->vert[0]->set(p0, 4*index + 0);
|
|
this->vert[1]->set(p1, 4*index + 1);
|
|
this->vert[2]->set(p2, 4*index + 2);
|
|
this->vert[3]->set(p3, 4*index + 3);
|
|
|
|
// determine the tetrahedron's type
|
|
double t = dot(p3 - p0, cross(p1 - p0, p2 - p0));
|
|
|
|
if (t < -1.0*GCTL_ZERO)
|
|
{
|
|
vec_order = regular_order;
|
|
return;
|
|
}
|
|
|
|
if (t > GCTL_ZERO)
|
|
{
|
|
vec_order = inverse_order;
|
|
return;
|
|
}
|
|
|
|
throw invalid_argument("invalid vertex positions. From type_tetrahedron::set(...)");
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::set(const point3ds &ps0, const point3ds &ps1,
|
|
const point3ds &ps2, const point3ds &ps3, int index)
|
|
{
|
|
if (index < 0)
|
|
{
|
|
throw out_of_range("Invalid index number, From type_tetrahedron::set(...)");
|
|
}
|
|
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
this->vert[i] = new vertex3dc;
|
|
}
|
|
this->self_host = true;
|
|
|
|
point3dc p0 = ps0.s2c();
|
|
point3dc p1 = ps1.s2c();
|
|
point3dc p2 = ps2.s2c();
|
|
point3dc p3 = ps3.s2c();
|
|
|
|
this->id = index;
|
|
this->vert[0]->set(p0, 4*index + 0);
|
|
this->vert[1]->set(p1, 4*index + 1);
|
|
this->vert[2]->set(p2, 4*index + 2);
|
|
this->vert[3]->set(p3, 4*index + 3);
|
|
|
|
// determine the tetrahedron's type
|
|
double t = dot(p3 - p0, cross(p1 - p0, p2 - p0));
|
|
|
|
if (t < -1.0*GCTL_ZERO)
|
|
{
|
|
vec_order = regular_order;
|
|
return;
|
|
}
|
|
|
|
if (t > GCTL_ZERO)
|
|
{
|
|
vec_order = inverse_order;
|
|
return;
|
|
}
|
|
|
|
throw invalid_argument("invalid vertex positions. From type_tetrahedron::set(...)");
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::set(vertex3dc &vert0, vertex3dc &vert1,
|
|
vertex3dc &vert2, vertex3dc &vert3, int index)
|
|
{
|
|
if (index < 0)
|
|
{
|
|
throw invalid_argument("Invalid index, From type_tetrahedron::set(...)");
|
|
}
|
|
|
|
this->id = index;
|
|
this->vert[0] = &vert0;
|
|
this->vert[1] = &vert1;
|
|
this->vert[2] = &vert2;
|
|
this->vert[3] = &vert3;
|
|
|
|
// determine the tetrahedron's type
|
|
double t = dot(vert3 - vert0, cross(vert1 - vert0, vert2 - vert0));
|
|
|
|
if (t < -1.0*GCTL_ZERO)
|
|
{
|
|
vec_order = regular_order;
|
|
return;
|
|
}
|
|
|
|
if (t > GCTL_ZERO)
|
|
{
|
|
vec_order = inverse_order;
|
|
return;
|
|
}
|
|
|
|
throw invalid_argument("invalid vertex positions. From type_tetrahedron::set(...)");
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::set(vertex3dc *vertp0, vertex3dc *vertp1,
|
|
vertex3dc *vertp2, vertex3dc *vertp3, int index)
|
|
{
|
|
if (index < 0)
|
|
{
|
|
throw invalid_argument("Invalid index, From type_tetrahedron::set(...)");
|
|
}
|
|
|
|
this->id = index;
|
|
this->vert[0] = vertp0;
|
|
this->vert[1] = vertp1;
|
|
this->vert[2] = vertp2;
|
|
this->vert[3] = vertp3;
|
|
|
|
// determine the tetrahedron's type
|
|
double t = dot(*vertp3 - *vertp0, cross(*vertp1 - *vertp0, *vertp2 - *vertp0));
|
|
|
|
if (t < -1.0*GCTL_ZERO)
|
|
{
|
|
vec_order = regular_order;
|
|
return;
|
|
}
|
|
|
|
if (t > GCTL_ZERO)
|
|
{
|
|
vec_order = inverse_order;
|
|
return;
|
|
}
|
|
|
|
throw invalid_argument("invalid vertex positions. From type_tetrahedron::set(...)");
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::reset()
|
|
{
|
|
entity<vertex3dc, 4, A>::reset();
|
|
vec_order = nullptr;
|
|
neigh[0] = neigh[1] = neigh[2] = neigh[3] = nullptr;
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::set_neighbor(type_tetrahedron<A> &nei0, type_tetrahedron<A> &nei1,
|
|
type_tetrahedron<A> &nei2, type_tetrahedron<A> &nei3)
|
|
{
|
|
if(!is_adjoined(nei0) || !is_adjoined(nei1) || !is_adjoined(nei2) || !is_adjoined(nei3))
|
|
{
|
|
throw invalid_argument("Invalid neighbors. From type_etrahedron::set_neighbor(...)");
|
|
}
|
|
|
|
neigh[0] = &nei0;
|
|
neigh[1] = &nei1;
|
|
neigh[2] = &nei2;
|
|
neigh[3] = &nei3;
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
vertex3dc *type_tetrahedron<A>::get(unsigned int i, unsigned int j) const
|
|
{
|
|
if (i > 3 || j > 2)
|
|
{
|
|
throw out_of_range("Invalid facet or vertex index. From type_tetrahedron::get(...)");
|
|
}
|
|
|
|
if (vec_order == nullptr)
|
|
{
|
|
throw domain_error("object not initialized. From type_tetrahedron::get(...)");
|
|
}
|
|
|
|
if (this->vert[0] == nullptr || this->vert[1] == nullptr ||
|
|
this->vert[2] == nullptr || this->vert[3] == nullptr)
|
|
{
|
|
throw domain_error("Invalid pointer. From type_tetrahedron::get(...)");
|
|
}
|
|
|
|
return this->vert[vec_order[3*i+j]];
|
|
}
|
|
|
|
template <typename A>
|
|
vertex3dc *type_tetrahedron<A>::fget(unsigned int i, unsigned int j) const
|
|
{
|
|
return this->vert[vec_order[3*i+j]];
|
|
}
|
|
|
|
template <typename A>
|
|
void type_tetrahedron<A>::deter_vert_order()
|
|
{
|
|
if (this->vert[0] == nullptr || this->vert[1] == nullptr ||
|
|
this->vert[2] == nullptr || this->vert[3] == nullptr)
|
|
{
|
|
throw domain_error("Invalid pointer. From type_tetrahedron::deter_vert_order(...)");
|
|
}
|
|
|
|
// determine the tetrahedron's type
|
|
double t = dot(*this->vert[3] - *this->vert[0],
|
|
cross(*this->vert[1] - *this->vert[0], *this->vert[2] - *this->vert[0]));
|
|
|
|
if (t < -1.0*GCTL_ZERO)
|
|
{
|
|
vec_order = regular_order;
|
|
return;
|
|
}
|
|
|
|
if (t > GCTL_ZERO)
|
|
{
|
|
vec_order = inverse_order;
|
|
return;
|
|
}
|
|
|
|
throw invalid_argument("invalid vertex positions. From type_tetrahedron::deter_vert_order(...)");
|
|
return;
|
|
}
|
|
|
|
template <typename A>
|
|
bool type_tetrahedron<A>::is_adjoined(type_tetrahedron<A> &nei)
|
|
{
|
|
int joined_count = 0;
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
for (int j = 0; j < 4; ++j)
|
|
{
|
|
if (this->vert[i] == nei.vert[j])
|
|
{
|
|
joined_count++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (joined_count == 3) return true;
|
|
return false;
|
|
}
|
|
|
|
template <typename A>
|
|
double type_tetrahedron<A>::volume()
|
|
{
|
|
// calculate element's volume
|
|
point3dc va, vb, vc;
|
|
va = *this->vert[1] - *this->vert[0];
|
|
vb = *this->vert[2] - *this->vert[0];
|
|
vc = *this->vert[3] - *this->vert[0];
|
|
return GCTL_FABS(dot(cross(va,vb),vc))/6.0;
|
|
}
|
|
|
|
template <typename A>
|
|
point3dc type_tetrahedron<A>::center()
|
|
{
|
|
point3dc vc;
|
|
vc = 0.25*(*this->vert[0] + *this->vert[1] + *this->vert[2] + *this->vert[3]);
|
|
return vc;
|
|
}
|
|
|
|
template <typename A>
|
|
void copy_type_tetrahedron(type_tetrahedron<A> *tar, const type_tetrahedron<A> *src)
|
|
{
|
|
copy_entity(tar, src);
|
|
|
|
tar->vec_order = src->vec_order;
|
|
|
|
for (size_t i = 0; i < 4; i++)
|
|
{
|
|
tar->neigh[i] = src->neigh[i];
|
|
}
|
|
return;
|
|
}
|
|
|
|
template <typename A, typename B>
|
|
void copy_type_tetrahedron(type_tetrahedron<A> *tar, const type_tetrahedron<B> *src)
|
|
{
|
|
copy_entity(tar, src);
|
|
|
|
tar->vec_order = src->vec_order;
|
|
return;
|
|
}
|
|
}
|
|
|
|
#endif // _GCTL_TETRAHEDRON_H
|