/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * 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 . * * 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_SHAPEFUNC_H #define _GCTL_SHAPEFUNC_H #include "../gctl_config.h" #include "../core/enum.h" #include "../core/exceptions.h" #include "../maths/mathfunc.h" namespace gctl { enum eshape_value_e { VecValue, // Original vector components VecGradient, // Gradients of the vector components VecCurl, // Curl components }; /** * @brief The 1st order (linear) shape functions */ class linear_sf { public: linear_sf(){} virtual ~linear_sf(){} /** * @brief Shape function defined on a 1D line. * * This function has two nodes on the left and right ends, respectively. * index for the left node function is 0. And the node's index is 1. * 0-----------1-----> x * * @param[in] x independent variable * @param[in] x1 low bound of the shape function * @param[in] x2 high bound of the shape function * @param[in] idx Node index * @param[in] vt Calculation target. Could be Value or Dx * * @return Value of the shape function */ double line(double x, double x1, double x2, size_t idx, gradient_type_e vt = Value); /** * @brief Shape function defined on a 2D rectangle. * * This function has four nodes on the corners, respectively. * index of the nodes are as show * * y * | * | * 3----------2 * | | * | | * | | * | | * 0----------1--->x * * @param[in] x independent variable x * @param[in] y independent variable y * @param[in] x1 low bound of x * @param[in] x2 high bound of x * @param[in] y1 low bound of y * @param[in] y2 high bound of y * @param[in] idx Node index * @param[in] vt Calculation target. Could be Value, Dx or Dy * * @return Value of the shape function */ double quad(double x, double y, double x1, double x2, double y1, double y2, size_t idx, gradient_type_e vt = Value); /** * @brief Shape function defined on a 2D triangle. * * This function has three nodes on the corners, respectively. * index of the nodes are as show * * y * | * | * 2 * | \ * | \ * | \ * | \ * | \ * | \ * 0-------------1--->x * * @param[in] x independent variable x * @param[in] y independent variable y * @param[in] x1 x coordinate of the first node of the triangle * @param[in] x2 x coordinate of the second node of the triangle * @param[in] x3 x coordinate of the third node of the triangle * @param[in] y1 y coordinate of the first node of the triangle * @param[in] y2 y coordinate of the second node of the triangle * @param[in] y3 y coordinate of the third node of the triangle * @param[in] idx Node index * @param[in] vt Calculation target. Could be Value, Dx or Dy * * @return Value of the shape function */ double triangle(double x, double y, double x1, double x2, double x3, double y1, double y2, double y3, size_t idx, gradient_type_e vt = Value); /** * @brief Interpolate using the triangular linear shape function * * This function has three nodes on the corners, respectively. * index of the nodes are as show * * y * | * | * 2 * | \ * | \ * | \ * | \ * | \ * | \ * 0-------------1--->x * * @param[in] x independent variable x * @param[in] y independent variable y * @param[in] x1 x coordinate of the first node of the triangle * @param[in] x2 x coordinate of the second node of the triangle * @param[in] x3 x coordinate of the third node of the triangle * @param[in] y1 y coordinate of the first node of the triangle * @param[in] y2 y coordinate of the second node of the triangle * @param[in] y3 y coordinate of the third node of the triangle * @param[in] v1 value at the first node of the triangle * @param[in] v2 value at the second node of the triangle * @param[in] v3 value at of the third node of the triangle * * @return Value of the interpolation */ double triangle_interpolate(double x, double y, double x1, double x2, double x3, double y1, double y2, double y3, double v1, double v2, double v3, gradient_type_e vt); /** * @brief Shape function defined on a 3D cube. * * This function has eight nodes on the corners, respectively. * index of the nodes are as show * z * | * | 7----------6 * | /| /| * | / | / | * |/ | / | * 4----------5 | * | 3 | 2 * | /-------|--/ * | / | / * 0----------1--------> x * * @param[in] x independent variable x * @param[in] y independent variable y * @param[in] z independent variable z * @param[in] x1 low bound of x * @param[in] x2 high bound of x * @param[in] y1 low bound of y * @param[in] y2 high bound of y * @param[in] z1 low bound of z * @param[in] z2 high bound of z * @param[in] idx Node index * @param[in] vt Calculation target. Could be Value, Dx, Dy or Dz * * @return Value of the shape function */ double cube(double x, double y, double z, double x1, double x2, double y1, double y2, double z1, double z2, size_t idx, gradient_type_e vt = Value); /** * @brief Shape function defined on a 3D tetrahedron. * * This function has four nodes on the corners, respectively. * index of the nodes are as show * zta * | * 3 * |\ * | \ * | \ * | \ * | \ * | \ * 0------------2---> eta * | / * | / * | / * | / * | / * | / * |/ * 1 * | * ksi * * @param[in] ksi local coordinate ksi * @param[in] eta local coordinate eta * @param[in] zta local coordinate zta * @param[in] x1 x coordinate of the first node of the tetrahedron * @param[in] x2 x coordinate of the second node of the tetrahedron * @param[in] x3 x coordinate of the third node of the tetrahedron * @param[in] x4 x coordinate of the fourth node of the tetrahedron * @param[in] y1 y coordinate of the first node of the tetrahedron * @param[in] y2 y coordinate of the second node of the tetrahedron * @param[in] y3 y coordinate of the third node of the tetrahedron * @param[in] y4 y coordinate of the fourth node of the tetrahedron * @param[in] z1 z coordinate of the first node of the tetrahedron * @param[in] z2 z coordinate of the second node of the tetrahedron * @param[in] z3 z coordinate of the third node of the tetrahedron * @param[in] z4 z coordinate of the fourth node of the tetrahedron * @param[in] idx Node index * @param[in] vt Calculation target. Could be Value, Dx, Dy or Dz * * @return Value of the shape function */ double tetrahedron(double ksi, double eta, double zta, double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4, size_t idx, gradient_type_e vt = Value); /** * @brief transform local coordinates to global ones within a tetrahedron element * * @param[out] x global coordinate x * @param[out] y global coordinate y * @param[out] z global coordinate z * @param[in] ksi local coordinate ksi * @param[in] eta local coordinate eta * @param[in] zta local coordinate zta * @param[in] x1 x coordinate of the first node of the tetrahedron * @param[in] x2 x coordinate of the second node of the tetrahedron * @param[in] x3 x coordinate of the third node of the tetrahedron * @param[in] x4 x coordinate of the fourth node of the tetrahedron * @param[in] y1 y coordinate of the first node of the tetrahedron * @param[in] y2 y coordinate of the second node of the tetrahedron * @param[in] y3 y coordinate of the third node of the tetrahedron * @param[in] y4 y coordinate of the fourth node of the tetrahedron * @param[in] z1 z coordinate of the first node of the tetrahedron * @param[in] z2 z coordinate of the second node of the tetrahedron * @param[in] z3 z coordinate of the third node of the tetrahedron * @param[in] z4 z coordinate of the fourth node of the tetrahedron */ void local2global_tetrahedron(double &x, double &y, double &z, double ksi, double eta, double zta, double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4); /** * @brief transform global coordinates to local ones within a tetrahedron element * * @param[out] ksi local coordinate ksi * @param[out] eta local coordinate eta * @param[out] zta local coordinate zta * @param[in] x global coordinate x * @param[in] y global coordinate y * @param[in] z global coordinate z * @param[in] x1 x coordinate of the first node of the tetrahedron * @param[in] x2 x coordinate of the second node of the tetrahedron * @param[in] x3 x coordinate of the third node of the tetrahedron * @param[in] x4 x coordinate of the fourth node of the tetrahedron * @param[in] y1 y coordinate of the first node of the tetrahedron * @param[in] y2 y coordinate of the second node of the tetrahedron * @param[in] y3 y coordinate of the third node of the tetrahedron * @param[in] y4 y coordinate of the fourth node of the tetrahedron * @param[in] z1 z coordinate of the first node of the tetrahedron * @param[in] z2 z coordinate of the second node of the tetrahedron * @param[in] z3 z coordinate of the third node of the tetrahedron * @param[in] z4 z coordinate of the fourth node of the tetrahedron */ void global2local_tetrahedron(double &ksi, double &eta, double &zta, double x, double y, double z, double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4); /** * @brief Interpolate using the tetrahedral linear shape function * * @param[in] x global coordinate x * @param[in] y global coordinate y * @param[in] z global coordinate z * @param[in] x1 x coordinate of the first node of the tetrahedron * @param[in] x2 x coordinate of the second node of the tetrahedron * @param[in] x3 x coordinate of the third node of the tetrahedron * @param[in] x4 x coordinate of the fourth node of the tetrahedron * @param[in] y1 y coordinate of the first node of the tetrahedron * @param[in] y2 y coordinate of the second node of the tetrahedron * @param[in] y3 y coordinate of the third node of the tetrahedron * @param[in] y4 y coordinate of the fourth node of the tetrahedron * @param[in] z1 z coordinate of the first node of the tetrahedron * @param[in] z2 z coordinate of the second node of the tetrahedron * @param[in] z3 z coordinate of the third node of the tetrahedron * @param[in] z4 z coordinate of the fourth node of the tetrahedron * @param vt Value type */ double tetrahedron_interpolate(double x, double y, double z, double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4, double v1, double v2, double v3, double v4, gradient_type_e vt = Value); /** * @brief Shape function defined on a 3D triangular prism. * * This function has six nodes on the corners, respectively. Index of the nodes are as show. * * zta eta * | / 3(0,1,1) * | / /|\ * |/ / | \ * | / | \ * | / / \0(0,1,0) * | 4|---------| 5 * | | / \ | * | |/-------\| * | 1(0,0,0) 2(1,0,0) * O-----------------> ksi * * @param[in] ksi local coordinate ksi * @param[in] eta local coordinate eta * @param[in] zta local coordinate zta * @param[in] xs x coordinates of the nodes of the prism, must be in length of six * @param[in] ys y coordinates of the nodes of the prism, must be in length of six * @param[in] zs z coordinates of the nodes of the prism, must be in length of six * @param[in] idx node index of the calculation * @param[in] vt Calculation target. Could be Value, Dx, Dy or Dz * @return Value of the shape function */ double triprism(double ksi, double eta, double zta, double *xs, double *ys, double *zs, size_t idx, gradient_type_e vt = Value); /** * @brief transform local coordinates to global ones within a triangular prism element * * @param[out] x global coordinate x * @param[out] y global coordinate y * @param[out] z global coordinate z * @param[in] ksi local coordinate ksi * @param[in] eta local coordinate eta * @param[in] zta local coordinate zta * @param[in] xs x coordinates of the nodes of the prism, must be in length of six * @param[in] ys y coordinates of the nodes of the prism, must be in length of six * @param[in] zs z coordinates of the nodes of the prism, must be in length of six */ void local2global_triprism(double &x, double &y, double &z, double ksi, double eta, double zta, double *xs, double *ys, double *zs); }; class linear_esf { public: linear_esf(){} virtual ~linear_esf(){} /** * @brief Edge based vector shape function defined on a 2D triangle. * * This function has three nodes on the corners, respectively. * Indice of the nodes and edges are as show. * * y * | * | * 2 * | \ * | \ * | \ * | \ * | \ * | \ * 0-------------1--->x * edge0: 0->1 * edge1: 1->2 * edge2: 2->0 * * @param[in] x Evaluation position x * @param[in] y Evaluation position y * @param[in] x1 x coordinate of the first node of the triangle * @param[in] x2 x coordinate of the second node of the triangle * @param[in] x3 x coordinate of the third node of the triangle * @param[in] y1 y coordinate of the first node of the triangle * @param[in] y2 y coordinate of the second node of the triangle * @param[in] y3 y coordinate of the third node of the triangle * @param[in] idx Edge index * @param[in] vt Calculation target. Could be Value, Gradient or Curl * @param[in] ot Ordering direction of the edge's nodes * @param xval Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V) (Curl) * @param yval Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V) (Curl) */ void triangle(double x, double y, double x1, double x2, double x3, double y1, double y2, double y3, size_t idx, eshape_value_e vt, edge_orient_e ot, double &xval, double &yval); /** * @brief Edge based vector shape function defined on a 3D tetrahedron. * * This function has four nodes on the corners, respectively. * Indice of the nodes and edges are as show. * * q * | * 3 * |\ * | \ * | \ * | \ * | \ * | \ * 0------------2---> v * | / * | / * | / * | / * | / * | / * |/ * 1 * | * u * * edge0: 0->1 * edge1: 0->2 * edge2: 0->3 * edge3: 1->2 * edge4: 1->3 (or 3->1 in some other implementations, Does not really matter ^,^) * edge5: 2->3 * * @param[in] x Evaluation position x * @param[in] y Evaluation position y * @param[in] z Evaluation position z * @param[in] x1 x coordinate of the first node of the tetrahedron * @param[in] x2 x coordinate of the second node of the tetrahedron * @param[in] x3 x coordinate of the third node of the tetrahedron * @param[in] x4 x coordinate of the fourth node of the tetrahedron * @param[in] y1 y coordinate of the first node of the tetrahedron * @param[in] y2 y coordinate of the second node of the tetrahedron * @param[in] y3 y coordinate of the third node of the tetrahedron * @param[in] y4 y coordinate of the fourth node of the tetrahedron * @param[in] z1 z coordinate of the first node of the tetrahedron * @param[in] z2 z coordinate of the second node of the tetrahedron * @param[in] z3 z coordinate of the third node of the tetrahedron * @param[in] z4 z coordinate of the fourth node of the tetrahedron * @param[in] idx Edge index * @param[in] vt Calculation target. Could be Value, Curl or Div * @param[in] ot Ordering direction of the edge's nodes * @param xval Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V)_x (Curl) * @param yval Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V)_y (Curl) * @param zval Evaluated shape function. Equals to Vz (Value), dVz/dx (Gradient) or curl(V)_z (Curl) * @param xval2 Evaluated shape function. Equals to dVx/dz (Gradient) * @param yval2 Evaluated shape function. Equals to dVy/dz (Gradient) * @param zval2 Evaluated shape function. Equals to dVz/dy (Gradient) */ void tetrahedron(double x, double y, double z, double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4, size_t idx, eshape_value_e vt, edge_orient_e ot, double &xval, double &yval, double &zval, double &xval2, double &yval2, double &zval2); /** * @brief Edge based vector shape function defined on a 3D cube. * * This function has eight nodes on the corners, respectively. * index of the nodes and edges are as show. * z * | * | 7----------6 * | /| /| * | / | / | * |/ | / | * 4----------5 | * | 3 -----|-- 2 * | / | / * | / | / * 0----------1--------> x * * edge0: 0->1 * edge1: 3->2 * edge2: 4->5 * edge3: 7->6 * edge4: 0->3 * edge5: 1->2 * edge6: 4->7 * edge7: 5->6 * edge8: 0->4 * edge9: 1->5 * edge10: 3->7 * edge11: 2->6 * * @param[in] x Evaluation variable x * @param[in] y Evaluation variable y * @param[in] z Evaluation variable z * @param[in] x1 low bound of x * @param[in] x2 high bound of x * @param[in] y1 low bound of y * @param[in] y2 high bound of y * @param[in] z1 low bound of z * @param[in] z2 high bound of z * @param[in] idx Edge index * @param[in] vt Calculation target. Could be Value, Curl or Div * @param[in] ot Ordering direction of the edge's nodes * @param xval Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V)_x (Curl) * @param yval Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V)_y (Curl) * @param zval Evaluated shape function. Equals to Vz (Value), dVz/dx (Gradient) or curl(V)_z (Curl) * @param xval2 Evaluated shape function. Equals to dVx/dz (Gradient) * @param yval2 Evaluated shape function. Equals to dVy/dz (Gradient) * @param zval2 Evaluated shape function. Equals to dVz/dy (Gradient) */ void cube(double x, double y, double z, double x1, double x2, double y1, double y2, double z1, double z2, size_t idx, eshape_value_e vt, edge_orient_e ot, double &xval, double &yval, double &zval, double &xval2, double &yval2, double &zval2); }; } #endif // _GCTL_SHAPEFUNC_H