tetra_element_linear.h

Go to the documentation of this file.

// 
//  This file is part of the FFEA simulation package
//  
//  Copyright (c) by the Theory and Development FFEA teams,
//  as they appear in the README.md file. 
// 
//  FFEA 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 3 of the License, or
//  (at your option) any later version.
// 
//  FFEA is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
// 
//  You should have received a copy of the GNU General Public License
//  along with FFEA.  If not, see <http://www.gnu.org/licenses/>.
// 
//  To help us fund FFEA development, we humbly ask that you cite 
//  the research papers on the package.
//

/*
 * tetra_element_linear.hpp
 */

#ifndef TETRA_ELEMENT_LINEAR_H_INCLUDED
#define TETRA_ELEMENT_LINEAR_H_INCLUDED

#include <stddef.h>
#include <stdio.h>
#include <map>
#include "mat_vec_fns.h"
#include "SimulationParams.h"
#include "FFEA_return_codes.h"
#include "RngStream.h"
#include "SecondOrderFunctions.h"
#include "PoissonMatrixQuadratic.h"
#include "MassMatrixLinear.h"
#include "MassMatrixQuadratic.h"

class Blob;

#define NUM_NODES_LINEAR_TET 4
#define NUM_NODES_QUADRATIC_TET 10

/*
 * A preprocessor macro generating code that calculates a submatrix of the bulk viscous
 * matrix. 'I' specifies what position in the matrix 'V' this 4x4 submatrix should be
 * inserted to (position I,I in the matrix).
 *
 * N.B. This macro only deals with submatrices *on the diagonal* of the bulk viscous matrix.
 */
#define BULK_VISCOUS_SUBMATRIX_DIAG(V,I) \
    V[I + 0][I + 0] = dpsi[I + 0] * dpsi[I + 0] * (A + B); \
    V[I + 1][I + 1] = dpsi[I + 1] * dpsi[I + 1] * (A + B); \
    V[I + 2][I + 2] = dpsi[I + 2] * dpsi[I + 2] * (A + B); \
    V[I + 3][I + 3] = dpsi[I + 3] * dpsi[I + 3] * (A + B); \
    V[I + 0][I + 1] = dpsi[I + 1] * dpsi[I + 0] * (A + B); \
    V[I + 0][I + 2] = dpsi[I + 2] * dpsi[I + 0] * (A + B); \
    V[I + 0][I + 3] = dpsi[I + 3] * dpsi[I + 0] * (A + B); \
    V[I + 1][I + 2] = dpsi[I + 1] * dpsi[I + 2] * (A + B); \
    V[I + 1][I + 3] = dpsi[I + 1] * dpsi[I + 3] * (A + B); \
    V[I + 2][I + 3] = dpsi[I + 2] * dpsi[I + 3] * (A + B);

/*
 * A preprocessor macro generating code that calculates a submatrix of the bulk viscous
 * matrix. 'I' and 'J' specify what position in the matrix 'V' this 4x4 submatrix should
 * be inserted to (position I,J in the matrix).
 *
 * N.B. This code will calculate submatrices anywhere in the bulk viscous matrix,
 * but the BULK_VISCOUS_SUBMATRIX_DIAG macro above provides simpler more efficient
 * code for the case in which the submatrix lies on the diagonal of the bulk viscous
 * matrix.
 */
#define BULK_VISCOUS_SUBMATRIX_OFFDIAG(V, I, J) \
    K[0][0] = dpsi[I] * dpsi[J]; \
    K[0][1] = dpsi[J + 1] * dpsi[I]; \
    K[0][2] = dpsi[J + 2] * dpsi[I]; \
    K[0][3] = dpsi[J + 3] * dpsi[I]; \
    \
    K[1][0] = dpsi[I + 1] * dpsi[J]; \
    K[1][1] = dpsi[I + 1] * dpsi[J + 1]; \
    K[1][2] = dpsi[I + 1] * dpsi[J + 2]; \
    K[1][3] = dpsi[I + 1] * dpsi[J + 3]; \
    \
    K[2][0] = dpsi[I + 2] * dpsi[J]; \
    K[2][1] = dpsi[I + 2] * dpsi[J + 1]; \
    K[2][2] = dpsi[I + 2] * dpsi[J + 2]; \
    K[2][3] = dpsi[I + 2] * dpsi[J + 3]; \
    \
    K[3][0] = dpsi[I + 3] * dpsi[J]; \
    K[3][1] = dpsi[I + 3] * dpsi[J + 1]; \
    K[3][2] = dpsi[I + 3] * dpsi[J + 2]; \
    K[3][3] = dpsi[I + 3] * dpsi[J + 3]; \
    \
    V[I + 0][J + 0] = (B + A) * K[0][0]; \
    V[I + 0][J + 1] = B * K[0][1] + A * K[1][0]; \
    V[I + 0][J + 2] = B * K[0][2] + A * K[2][0]; \
    V[I + 0][J + 3] = B * K[0][3] + A * K[3][0]; \
    \
    V[I + 1][J + 0] = B * K[1][0] + A * K[0][1]; \
    V[I + 1][J + 1] = (B + A) * K[1][1]; \
    V[I + 1][J + 2] = B * K[1][2] + A * K[2][1]; \
    V[I + 1][J + 3] = B * K[1][3] + A * K[3][1]; \
    \
    V[I + 2][J + 0] = B * K[2][0] + A * K[0][2]; \
    V[I + 2][J + 1] = B * K[2][1] + A * K[1][2]; \
    V[I + 2][J + 2] = (B + A) * K[2][2]; \
    V[I + 2][J + 3] = B * K[2][3] + A * K[3][2]; \
    \
    V[I + 3][J + 0] = B * K[3][0] + A * K[0][3]; \
    V[I + 3][J + 1] = B * K[3][1] + A * K[1][3]; \
    V[I + 3][J + 2] = B * K[3][2] + A * K[2][3]; \
    V[I + 3][J + 3] = (B + A) * K[3][3]; \

#define RAND(A, B) ((A) + ((B)-(A))*(rng[thread_id].RandU01()))

#define DPSI1_DX 0
#define DPSI2_DX 1
#define DPSI3_DX 2
#define DPSI4_DX 3
#define DPSI1_DY 4
#define DPSI2_DY 5
#define DPSI3_DY 6
#define DPSI4_DY 7
#define DPSI1_DZ 8
#define DPSI2_DZ 9
#define DPSI3_DZ 10
#define DPSI4_DZ 11

class tetra_element_linear {
public:

    tetra_element_linear();

    /* Properties of this element */
    scalar rho; 
    scalar A; 
    scalar B; 
    scalar G; 
    scalar E; 
    scalar dielectric;
    scalar mass;

    mesh_node *n[NUM_NODES_QUADRATIC_TET];

    vector12 dpsi;

    upper_triangular_matrix4 del2;

    PoissonMatrixQuadratic K_alpha;

    vector3 node_force[NUM_NODES_QUADRATIC_TET];

    scalar vol_0;

    scalar vol;

    matrix3 F_ij;

    scalar internal_stress_mag;

    matrix3 J_inv_0;

    matrix12 viscosity_matrix;

    /* Blob this element is a part of */
    Blob *daddy_blob;

    int index;

    vector3 centroid;

    scalar * get_K_alpha_element_mem_loc(int ni, int nj);

    void calculate_K_alpha();

    void construct_element_mass_matrix(MassMatrixQuadratic *M_alpha);
    void construct_element_mass_matrix(MassMatrixLinear *M_alpha);

    void add_K_alpha(scalar *K, int num_nodes);

    void get_grad_phi_at_stu(vector3 &grad_phi, scalar s, scalar t, scalar u);

    void calculate_electrostatic_forces();

    void calculate_jacobian(matrix3 J);

    void calc_deformation(matrix3 J);

    void calc_elastic_force_vector(vector12 F);

    int calc_shape_function_derivatives_and_volume(matrix3 J);

    void print_structural_details();

    scalar calc_volume();

    void create_viscosity_matrix();

    /*
     *
     */
    void add_shear_elastic_stress(matrix3 J, matrix3 stress);

    /*
     *
     */
    void add_bulk_elastic_stress(matrix3 stress);
    void add_bulk_elastic_stress_OLD(matrix3 stress);

    void add_fluctuating_stress(SimulationParams *params, RngStream rng[], matrix3 stress, int thread_id);

    void apply_stress_tensor(matrix3 stress, vector12 du);

    void get_element_velocity_vector(vector12 v);

    void add_element_force_vector(vector12 force);

    void add_force_to_node(int i, vector3 *f);

    int what_node_is_this(int index);

    void print();
    void print_viscosity_matrix();

    void apply_element_mass_matrix(vector12 du);

    void volume_coord_to_xyz(scalar eta0, scalar eta1, scalar eta2, scalar eta3, vector3 *r);

    void zero_force();

    void linearise_element();

    void calc_centroid();

    int get_opposite_node(int n1, int n2, int n3); 

    scalar length_of_longest_edge();

private:

    scalar last_det;

    void calc_del2_matrix();

    void add_diffusion_matrix(matrix12 V);

    struct tetrahedron_gauss_point {
        scalar W;
        scalar eta[4];
    };

};

#endif