Homec4science
Diffusion lammps (master) Browse

lammps/lib/atc/FE_Quadrature.h
0966e14e73a7master

Last Change
Raw File

FE_Quadrature.h
View Options

#ifndef FE_QUADRATURE_H
#define FE_QUADRATURE_H
#include <vector>
#include "MatrixLibrary.h"
#include "ATC_TypeDefs.h"
namespace ATC {
/**
* @class FE_Quadrature
* @brief Stores quadrature schemes
*/
struct FE_Quadrature {
int numIPs;
int numFaceIPs;
DENS_MAT ipCoords;
std::vector<DENS_MAT> ipFaceCoords;
DENS_MAT ipFace2DCoords;
DENS_VEC ipWeights;
DENS_VEC ipFaceWeights;
FE_Quadrature(FeEltGeometry geo, FeIntQuadrature quad)
{
if (geo == HEXA) {
switch (quad) {
// Degenerate 1 point quadrature
case GAUSS1: {
// Set number of IPs and face IPs
numIPs = 1;
numFaceIPs = 1;
// Size matrices and vectors thereof accordingly
ipCoords.resize(3,numIPs);
ipFaceCoords.assign(6,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(2,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
// "Matrix" of integration point location
ipCoords(0,0) = 0.0;
ipCoords(1,0) = 0.0;
ipCoords(2,0) = 0.0;
// Integration point for each face
ipFaceCoords[0](0,0) = -1.0; ipFaceCoords[3](0,0) = 0.0;
ipFaceCoords[0](1,0) = 0.0; ipFaceCoords[3](1,0) = 1.0;
ipFaceCoords[0](2,0) = 0.0; ipFaceCoords[3](2,0) = 0.0;
//
ipFaceCoords[1](0,0) = 1.0; ipFaceCoords[4](0,0) = 0.0;
ipFaceCoords[1](1,0) = 0.0; ipFaceCoords[4](1,0) = 0.0;
ipFaceCoords[1](2,0) = 0.0; ipFaceCoords[4](2,0) = -1.0;
//
ipFaceCoords[2](0,0) = 0.0; ipFaceCoords[5](0,0) = 0.0;
ipFaceCoords[2](1,0) = -1.0; ipFaceCoords[5](1,0) = 0.0;
ipFaceCoords[2](2,0) = 0.0; ipFaceCoords[5](2,0) = 1.0;
// 2D integration scheme for the faces
ipFace2DCoords(0,0) = 0.0;
ipFace2DCoords(1,0) = 0.0;
// Integration point weights
ipWeights = 8.0;
// Face integration point weights
ipFaceWeights = 4.0;
break;
}
// 8 point quadratures
case NODAL:
case GAUSS2: {
// Set number of IPs and face IPs
numIPs = 8;
numFaceIPs = 4;
// Size matrices and vectors thereof accordingly
ipCoords.resize(3,numIPs);
ipFaceCoords.assign(6,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(2,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
// Dictates difference in node locations for nodal/GAUSS2
double a = 1.0/sqrt(3.0);
if (quad == NODAL) a = 1.0;
// Matrix of integration point locations & follows local
// conn
ipCoords(0,0) = -a; ipCoords(0,4) = -a;
ipCoords(1,0) = -a; ipCoords(1,4) = -a;
ipCoords(2,0) = -a; ipCoords(2,4) = a;
//
ipCoords(0,1) = a; ipCoords(0,5) = a;
ipCoords(1,1) = -a; ipCoords(1,5) = -a;
ipCoords(2,1) = -a; ipCoords(2,5) = a;
//
ipCoords(0,2) = a; ipCoords(0,6) = a;
ipCoords(1,2) = a; ipCoords(1,6) = a;
ipCoords(2,2) = -a; ipCoords(2,6) = a;
//
ipCoords(0,3) = -a; ipCoords(0,7) = -a;
ipCoords(1,3) = a; ipCoords(1,7) = a;
ipCoords(2,3) = -a; ipCoords(2,7) = a;
// Integration points by face
ipFaceCoords[0](0,0) = -1; ipFaceCoords[3](0,0) = -a;
ipFaceCoords[0](1,0) = -a; ipFaceCoords[3](1,0) = 1;
ipFaceCoords[0](2,0) = -a; ipFaceCoords[3](2,0) = -a;
//
ipFaceCoords[0](0,1) = -1; ipFaceCoords[3](0,1) = a;
ipFaceCoords[0](1,1) = a; ipFaceCoords[3](1,1) = 1;
ipFaceCoords[0](2,1) = -a; ipFaceCoords[3](2,1) = -a;
//
ipFaceCoords[0](0,2) = -1; ipFaceCoords[3](0,2) = a;
ipFaceCoords[0](1,2) = a; ipFaceCoords[3](1,2) = 1;
ipFaceCoords[0](2,2) = a; ipFaceCoords[3](2,2) = a;
//
ipFaceCoords[0](0,3) = -1; ipFaceCoords[3](0,3) = -a;
ipFaceCoords[0](1,3) = -a; ipFaceCoords[3](1,3) = 1;
ipFaceCoords[0](2,3) = a; ipFaceCoords[3](2,3) = a;
ipFaceCoords[1](0,0) = 1; ipFaceCoords[4](0,0) = -a;
ipFaceCoords[1](1,0) = -a; ipFaceCoords[4](1,0) = -a;
ipFaceCoords[1](2,0) = -a; ipFaceCoords[4](2,0) = -1;
//
ipFaceCoords[1](0,1) = 1; ipFaceCoords[4](0,1) = a;
ipFaceCoords[1](1,1) = a; ipFaceCoords[4](1,1) = -a;
ipFaceCoords[1](2,1) = -a; ipFaceCoords[4](2,1) = -1;
//
ipFaceCoords[1](0,2) = 1; ipFaceCoords[4](0,2) = a;
ipFaceCoords[1](1,2) = a; ipFaceCoords[4](1,2) = a;
ipFaceCoords[1](2,2) = a; ipFaceCoords[4](2,2) = -1;
//
ipFaceCoords[1](0,3) = 1; ipFaceCoords[4](0,3) = -a;
ipFaceCoords[1](1,3) = -a; ipFaceCoords[4](1,3) = a;
ipFaceCoords[1](2,3) = a; ipFaceCoords[4](2,3) = -1;
ipFaceCoords[2](0,0) = -a; ipFaceCoords[5](0,0) = -a;
ipFaceCoords[2](1,0) = -1; ipFaceCoords[5](1,0) = -a;
ipFaceCoords[2](2,0) = -a; ipFaceCoords[5](2,0) = 1;
//
ipFaceCoords[2](0,1) = a; ipFaceCoords[5](0,1) = a;
ipFaceCoords[2](1,1) = -1; ipFaceCoords[5](1,1) = -a;
ipFaceCoords[2](2,1) = -a; ipFaceCoords[5](2,1) = 1;
//
ipFaceCoords[2](0,2) = a; ipFaceCoords[5](0,2) = a;
ipFaceCoords[2](1,2) = -1; ipFaceCoords[5](1,2) = a;
ipFaceCoords[2](2,2) = a; ipFaceCoords[5](2,2) = 1;
//
ipFaceCoords[2](0,3) = -a; ipFaceCoords[5](0,3) = -a;
ipFaceCoords[2](1,3) = -1; ipFaceCoords[5](1,3) = a;
ipFaceCoords[2](2,3) = a; ipFaceCoords[5](2,3) = 1;
// Integration points for all faces ignoring the
// redundant dim
ipFace2DCoords(0,0) = -a; ipFace2DCoords(0,2) = a;
ipFace2DCoords(1,0) = -a; ipFace2DCoords(1,2) = a;
//
ipFace2DCoords(0,1) = a; ipFace2DCoords(0,3) = -a;
ipFace2DCoords(1,1) = -a; ipFace2DCoords(1,3) = a;
// Integration point weights
ipWeights = 1.0;
// Face integration point weights
ipFaceWeights = 1.0;
break;
}
// 6 point "face" quadrature
case FACE: {
printf("using face quad!\n");
// Set number of IPs and face IPs
numIPs = 6;
numFaceIPs = 4;
// Size matrices and vectors thereof accordingly
ipCoords.resize(3,numIPs);
ipFaceCoords.assign(6,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(2,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
// Use GAUSS2 for faces for now...
double a = 1.0/sqrt(3.0);
// Matrix of integration point locations
ipCoords(0,0) = 1.0; ipCoords(0,3) = 0.0;
ipCoords(1,0) = 0.0; ipCoords(1,3) = -1.0;
ipCoords(2,0) = 0.0; ipCoords(2,3) = 0.0;
//
ipCoords(0,1) = -1.0; ipCoords(0,4) = 0.0;
ipCoords(1,1) = 0.0; ipCoords(1,4) = 0.0;
ipCoords(2,1) = 0.0; ipCoords(2,4) = 1.0;
//
ipCoords(0,2) = 0.0; ipCoords(0,5) = 0.0;
ipCoords(1,2) = 1.0; ipCoords(1,5) = 0.0;
ipCoords(2,2) = 0.0; ipCoords(2,5) = -1.0;
// Integration points by face
ipFaceCoords[0](0,0) = -1; ipFaceCoords[3](0,0) = -a;
ipFaceCoords[0](1,0) = -a; ipFaceCoords[3](1,0) = 1;
ipFaceCoords[0](2,0) = -a; ipFaceCoords[3](2,0) = -a;
//
ipFaceCoords[0](0,1) = -1; ipFaceCoords[3](0,1) = -a;
ipFaceCoords[0](1,1) = a; ipFaceCoords[3](1,1) = 1;
ipFaceCoords[0](2,1) = -a; ipFaceCoords[3](2,1) = a;
//
ipFaceCoords[0](0,2) = -1; ipFaceCoords[3](0,2) = a;
ipFaceCoords[0](1,2) = a; ipFaceCoords[3](1,2) = 1;
ipFaceCoords[0](2,2) = a; ipFaceCoords[3](2,2) = a;
//
ipFaceCoords[0](0,3) = -1; ipFaceCoords[3](0,3) = a;
ipFaceCoords[0](1,3) = -a; ipFaceCoords[3](1,3) = 1;
ipFaceCoords[0](2,3) = a; ipFaceCoords[3](2,3) = -a;
ipFaceCoords[1](0,0) = 1; ipFaceCoords[4](0,0) = -a;
ipFaceCoords[1](1,0) = -a; ipFaceCoords[4](1,0) = -a;
ipFaceCoords[1](2,0) = -a; ipFaceCoords[4](2,0) = -1;
//
ipFaceCoords[1](0,1) = 1; ipFaceCoords[4](0,1) = a;
ipFaceCoords[1](1,1) = a; ipFaceCoords[4](1,1) = -a;
ipFaceCoords[1](2,1) = -a; ipFaceCoords[4](2,1) = -1;
//
ipFaceCoords[1](0,2) = 1; ipFaceCoords[4](0,2) = a;
ipFaceCoords[1](1,2) = a; ipFaceCoords[4](1,2) = a;
ipFaceCoords[1](2,2) = a; ipFaceCoords[4](2,2) = -1;
//
ipFaceCoords[1](0,3) = 1; ipFaceCoords[4](0,3) = -a;
ipFaceCoords[1](1,3) = -a; ipFaceCoords[4](1,3) = a;
ipFaceCoords[1](2,3) = a; ipFaceCoords[4](2,3) = -1;
ipFaceCoords[2](0,0) = -a; ipFaceCoords[5](0,0) = -a;
ipFaceCoords[2](1,0) = -1; ipFaceCoords[5](1,0) = -a;
ipFaceCoords[2](2,0) = -a; ipFaceCoords[5](2,0) = 1;
//
ipFaceCoords[2](0,1) = -a; ipFaceCoords[5](0,1) = a;
ipFaceCoords[2](1,1) = -1; ipFaceCoords[5](1,1) = -a;
ipFaceCoords[2](2,1) = a; ipFaceCoords[5](2,1) = 1;
//
ipFaceCoords[2](0,2) = a; ipFaceCoords[5](0,2) = a;
ipFaceCoords[2](1,2) = -1; ipFaceCoords[5](1,2) = a;
ipFaceCoords[2](2,2) = a; ipFaceCoords[5](2,2) = 1;
//
ipFaceCoords[2](0,3) = a; ipFaceCoords[5](0,3) = -a;
ipFaceCoords[2](1,3) = -1; ipFaceCoords[5](1,3) = a;
ipFaceCoords[2](2,3) = -a; ipFaceCoords[5](2,3) = 1;
// Integration points for all faces ignoring the
// redundant dim
ipFace2DCoords(0,0) = -a; ipFace2DCoords(0,2) = a;
ipFace2DCoords(1,0) = -a; ipFace2DCoords(1,2) = a;
//
ipFace2DCoords(0,1) = a; ipFace2DCoords(0,3) = -a;
ipFace2DCoords(1,1) = -a; ipFace2DCoords(1,3) = a;
// Integration point weights
ipWeights = 4.0/3.0;
// Face integration point weights
ipFaceWeights = 1.0;
break;
}
// 27 point quadrature
case GAUSS3: {
// Set number of IPs and face IPs
numIPs = 27;
numFaceIPs = 9;
// Size matrices and vectors thereof accordingly
ipCoords.resize(3,numIPs);
ipFaceCoords.assign(6,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(2,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
double a = sqrt(3.0/5.0);
// Matrix of integration point locations & follows local
// conn
ipCoords(0,0) = -a; ipCoords(0,16) = -a;
ipCoords(1,0) = -a; ipCoords(1,16) = -a;
ipCoords(2,0) = -a; ipCoords(2,16) = 0;
//
ipCoords(0,1) = a; ipCoords(0,17) = -a;
ipCoords(1,1) = -a; ipCoords(1,17) = a;
ipCoords(2,1) = -a; ipCoords(2,17) = 0;
//
ipCoords(0,2) = a; ipCoords(0,18) = a;
ipCoords(1,2) = a; ipCoords(1,18) = -a;
ipCoords(2,2) = -a; ipCoords(2,18) = 0;
//
ipCoords(0,3) = -a; ipCoords(0,19) = a;
ipCoords(1,3) = a; ipCoords(1,19) = a;
ipCoords(2,3) = -a; ipCoords(2,19) = 0;
ipCoords(0,4) = -a; ipCoords(0,20) = 0;
ipCoords(1,4) = -a; ipCoords(1,20) = 0;
ipCoords(2,4) = a; ipCoords(2,20) = -a;
ipCoords(0,5) = a; ipCoords(0,21) = 0;
ipCoords(1,5) = -a; ipCoords(1,21) = 0;
ipCoords(2,5) = a; ipCoords(2,21) = a;
ipCoords(0,6) = a; ipCoords(0,22) = 0;
ipCoords(1,6) = a; ipCoords(1,22) = -a;
ipCoords(2,6) = a; ipCoords(2,22) = 0;
ipCoords(0,7) = -a; ipCoords(0,23) = 0;
ipCoords(1,7) = a; ipCoords(1,23) = a;
ipCoords(2,7) = a; ipCoords(2,23) = 0;
ipCoords(0,8) = 0; ipCoords(0,24) = -a;
ipCoords(1,8) = -a; ipCoords(1,24) = 0;
ipCoords(2,8) = -a; ipCoords(2,24) = 0;
ipCoords(0,9) = 0; ipCoords(0,25) = a;
ipCoords(1,9) = -a; ipCoords(1,25) = 0;
ipCoords(2,9) = a; ipCoords(2,25) = 0;
ipCoords(0,10) = 0; ipCoords(0,26) = 0;
ipCoords(1,10) = a; ipCoords(1,26) = 0;
ipCoords(2,10) = -a; ipCoords(2,26) = 0;
ipCoords(0,11) = 0;
ipCoords(1,11) = a;
ipCoords(2,11) = a;
ipCoords(0,12) = -a;
ipCoords(1,12) = 0;
ipCoords(2,12) = -a;
ipCoords(0,13) = -a;
ipCoords(1,13) = 0;
ipCoords(2,13) = a;
ipCoords(0,14) = a;
ipCoords(1,14) = 0;
ipCoords(2,14) = -a;
ipCoords(0,15) = a;
ipCoords(1,15) = 0;
ipCoords(2,15) = a;
// Integration points by face
ipFaceCoords[0](0,0) = -1; ipFaceCoords[0](0,5) = -1;
ipFaceCoords[0](1,0) = -a; ipFaceCoords[0](1,5) = 0;
ipFaceCoords[0](2,0) = -a; ipFaceCoords[0](2,5) = a;
//
ipFaceCoords[0](0,1) = -1; ipFaceCoords[0](0,6) = -1;
ipFaceCoords[0](1,1) = a; ipFaceCoords[0](1,6) = -a;
ipFaceCoords[0](2,1) = -a; ipFaceCoords[0](2,6) = 0;
//
ipFaceCoords[0](0,2) = -1; ipFaceCoords[0](0,7) = -1;
ipFaceCoords[0](1,2) = a; ipFaceCoords[0](1,7) = a;
ipFaceCoords[0](2,2) = a; ipFaceCoords[0](2,7) = 0;
//
ipFaceCoords[0](0,3) = -1; ipFaceCoords[0](0,8) = -1;
ipFaceCoords[0](1,3) = -a; ipFaceCoords[0](1,8) = 0;
ipFaceCoords[0](2,3) = a; ipFaceCoords[0](2,8) = 0;
//
ipFaceCoords[0](0,4) = -1;
ipFaceCoords[0](1,4) = 0;
ipFaceCoords[0](2,4) = -a;
ipFaceCoords[1](0,0) = 1; ipFaceCoords[1](0,5) = 1;
ipFaceCoords[1](1,0) = -a; ipFaceCoords[1](1,5) = 0;
ipFaceCoords[1](2,0) = -a; ipFaceCoords[1](2,5) = a;
//
ipFaceCoords[1](0,1) = 1; ipFaceCoords[1](0,6) = 1;
ipFaceCoords[1](1,1) = a; ipFaceCoords[1](1,6) = -a;
ipFaceCoords[1](2,1) = -a; ipFaceCoords[1](2,6) = 0;
//
ipFaceCoords[1](0,2) = 1; ipFaceCoords[1](0,7) = 1;
ipFaceCoords[1](1,2) = a; ipFaceCoords[1](1,7) = a;
ipFaceCoords[1](2,2) = a; ipFaceCoords[1](2,7) = 0;
//
ipFaceCoords[1](0,3) = 1; ipFaceCoords[1](0,8) = 1;
ipFaceCoords[1](1,3) = -a; ipFaceCoords[1](1,8) = 0;
ipFaceCoords[1](2,3) = a; ipFaceCoords[1](2,8) = 0;
//
ipFaceCoords[1](0,4) = 1;
ipFaceCoords[1](1,4) = 0;
ipFaceCoords[1](2,4) = -a;
ipFaceCoords[2](0,0) = -a; ipFaceCoords[2](0,5) = 0;
ipFaceCoords[2](1,0) = -1; ipFaceCoords[2](1,5) = -1;
ipFaceCoords[2](2,0) = -a; ipFaceCoords[2](2,5) = a;
//
ipFaceCoords[2](0,1) = -a; ipFaceCoords[2](0,6) = -a;
ipFaceCoords[2](1,1) = -1; ipFaceCoords[2](1,6) = -1;
ipFaceCoords[2](2,1) = a; ipFaceCoords[2](2,6) = 0;
//
ipFaceCoords[2](0,2) = a; ipFaceCoords[2](0,7) = a;
ipFaceCoords[2](1,2) = -1; ipFaceCoords[2](1,7) = -1;
ipFaceCoords[2](2,2) = a; ipFaceCoords[2](2,7) = 0;
//
ipFaceCoords[2](0,3) = a; ipFaceCoords[2](0,8) = 0;
ipFaceCoords[2](1,3) = -1; ipFaceCoords[2](1,8) = -1;
ipFaceCoords[2](2,3) = -a; ipFaceCoords[2](2,8) = 0;
//
ipFaceCoords[2](0,4) = 0;
ipFaceCoords[2](1,4) = -1;
ipFaceCoords[2](2,4) = -a;
ipFaceCoords[3](0,0) = -a; ipFaceCoords[3](0,5) = 0;
ipFaceCoords[3](1,0) = 1; ipFaceCoords[3](1,5) = 1;
ipFaceCoords[3](2,0) = -a; ipFaceCoords[3](2,5) = a;
//
ipFaceCoords[3](0,1) = -a; ipFaceCoords[3](0,6) = -a;
ipFaceCoords[3](1,1) = 1; ipFaceCoords[3](1,6) = 1;
ipFaceCoords[3](2,1) = a; ipFaceCoords[3](2,6) = 0;
//
ipFaceCoords[3](0,2) = a; ipFaceCoords[3](0,7) = a;
ipFaceCoords[3](1,2) = 1; ipFaceCoords[3](1,7) = 1;
ipFaceCoords[3](2,2) = a; ipFaceCoords[3](2,7) = 0;
//
ipFaceCoords[3](0,3) = a; ipFaceCoords[3](0,8) = 0;
ipFaceCoords[3](1,3) = 1; ipFaceCoords[3](1,8) = 1;
ipFaceCoords[3](2,3) = -a; ipFaceCoords[3](2,8) = 0;
//
ipFaceCoords[3](0,4) = 0;
ipFaceCoords[3](1,4) = 1;
ipFaceCoords[3](2,4) = -a;
ipFaceCoords[4](0,0) = -a; ipFaceCoords[4](0,5) = 0;
ipFaceCoords[4](1,0) = -a; ipFaceCoords[4](1,5) = a;
ipFaceCoords[4](2,0) = -1; ipFaceCoords[4](2,5) = -1;
//
ipFaceCoords[4](0,1) = a; ipFaceCoords[4](0,6) = -a;
ipFaceCoords[4](1,1) = -a; ipFaceCoords[4](1,6) = 0;
ipFaceCoords[4](2,1) = -1; ipFaceCoords[4](2,6) = -1;
//
ipFaceCoords[4](0,2) = a; ipFaceCoords[4](0,7) = a;
ipFaceCoords[4](1,2) = a; ipFaceCoords[4](1,7) = 0;
ipFaceCoords[4](2,2) = -1; ipFaceCoords[4](2,7) = -1;
//
ipFaceCoords[4](0,3) = -a; ipFaceCoords[4](0,8) = 0;
ipFaceCoords[4](1,3) = a; ipFaceCoords[4](1,8) = 0;
ipFaceCoords[4](2,3) = -1; ipFaceCoords[4](2,8) = -1;
//
ipFaceCoords[4](0,4) = 0;
ipFaceCoords[4](1,4) = -a;
ipFaceCoords[4](2,4) = -1;
ipFaceCoords[5](0,0) = -a; ipFaceCoords[5](0,5) = 0;
ipFaceCoords[5](1,0) = -a; ipFaceCoords[5](1,5) = a;
ipFaceCoords[5](2,0) = 1; ipFaceCoords[5](2,5) = 1;
//
ipFaceCoords[5](0,1) = a; ipFaceCoords[5](0,6) = -a;
ipFaceCoords[5](1,1) = -a; ipFaceCoords[5](1,6) = 0;
ipFaceCoords[5](2,1) = 1; ipFaceCoords[5](2,6) = 1;
//
ipFaceCoords[5](0,2) = a; ipFaceCoords[5](0,7) = a;
ipFaceCoords[5](1,2) = a; ipFaceCoords[5](1,7) = 0;
ipFaceCoords[5](2,2) = 1; ipFaceCoords[5](2,7) = 1;
//
ipFaceCoords[5](0,3) = -a; ipFaceCoords[5](0,8) = 0;
ipFaceCoords[5](1,3) = a; ipFaceCoords[5](1,8) = 0;
ipFaceCoords[5](2,3) = 1; ipFaceCoords[5](2,8) = 1;
//
ipFaceCoords[5](0,4) = 0;
ipFaceCoords[5](1,4) = -a;
ipFaceCoords[5](2,4) = 1;
// Integration points for all faces ignoring the
// redundant dim
ipFace2DCoords(0,0) = -a; ipFace2DCoords(0,5) = 0;
ipFace2DCoords(1,0) = -a; ipFace2DCoords(1,5) = a;
//
ipFace2DCoords(0,1) = a; ipFace2DCoords(0,6) = -a;
ipFace2DCoords(1,1) = -a; ipFace2DCoords(1,6) = 0;
//
ipFace2DCoords(0,2) = a; ipFace2DCoords(0,7) = a;
ipFace2DCoords(1,2) = a; ipFace2DCoords(1,7) = 0;
//
ipFace2DCoords(0,3) = -a; ipFace2DCoords(0,8) = 0;
ipFace2DCoords(1,3) = a; ipFace2DCoords(1,8) = 0;
//
ipFace2DCoords(0,4) = 0;
ipFace2DCoords(1,4) = -a;
// Integration point weights
for (int i=0; i<numIPs; ++i) {
if (i < 8) ipWeights[i] = 125.0/729.0;
else if (i < 20) ipWeights[i] = 200.0/729.0;
else if (i < 26) ipWeights[i] = 320.0/729.0;
else ipWeights[i] = 512.0/729.0;
}
// Face integration point weights
for (int i=0; i<numFaceIPs; ++i) {
if (i < 4) ipFaceWeights[i] = 25.0/81.0;
else if (i < 8) ipFaceWeights[i] = 40.0/81.0;
else ipFaceWeights[i] = 64.0/81.0;
}
break;
}
// Error
default: {
throw ATC_Error("Unrecognized quadrature type "
"for element type HEXA.");
}
}
} else if (geo == TETRA) {
switch (quad) {
// 4 point quadratures
case NODAL:
case GAUSS2: {
// Set number of IPs and face IPs
numIPs = 4;
numFaceIPs = 3;
// Size matrices and vectors thereof accordingly
ipCoords.resize(4,numIPs);
ipFaceCoords.assign(4,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(3,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
double v1, v2, a1, a2;
v1 = 0.585410196624969;
v2 = 0.138196601125011;
a1 = 2.0/3.0;
a2 = 1.0/6.0;
if (quad == NODAL) {
// Nodal quadrature
v1 = 1;
v2 = 0;
a1 = 1;
a2 = 0;
}
// Integration point coordinates
ipCoords.resize(4, numIPs);
ipCoords(0,0) = v2; ipCoords(0,2) = v2;
ipCoords(1,0) = v2; ipCoords(1,2) = v1;
ipCoords(2,0) = v2; ipCoords(2,2) = v2;
ipCoords(3,0) = v1; ipCoords(3,2) = v2;
//
ipCoords(0,1) = v1; ipCoords(0,3) = v2;
ipCoords(1,1) = v2; ipCoords(1,3) = v2;
ipCoords(2,1) = v2; ipCoords(2,3) = v1;
ipCoords(3,1) = v2; ipCoords(3,3) = v2;
// Integration points by face
// ...face 0 ...face 2
ipFaceCoords[0](0,0) = a1; ipFaceCoords[2](0,0) = a2;
ipFaceCoords[0](1,0) = a2; ipFaceCoords[2](1,0) = 0;
ipFaceCoords[0](2,0) = a2; ipFaceCoords[2](2,0) = a2;
//
ipFaceCoords[0](0,1) = a2; ipFaceCoords[2](0,1) = a1;
ipFaceCoords[0](1,1) = a1; ipFaceCoords[2](1,1) = 0;
ipFaceCoords[0](2,1) = a2; ipFaceCoords[2](2,1) = a2;
//
ipFaceCoords[0](0,2) = a2; ipFaceCoords[2](0,2) = a2;
ipFaceCoords[0](1,2) = a2; ipFaceCoords[2](1,2) = 0;
ipFaceCoords[0](2,2) = a1; ipFaceCoords[2](2,2) = a1;
// ...face 1 ...face 3
ipFaceCoords[1](0,0) = 0; ipFaceCoords[3](0,0) = a2;
ipFaceCoords[1](1,0) = a1; ipFaceCoords[3](1,0) = a2;
ipFaceCoords[1](2,0) = a2; ipFaceCoords[3](2,0) = 0;
//
ipFaceCoords[1](0,1) = 0; ipFaceCoords[3](0,1) = a2;
ipFaceCoords[1](1,1) = a2; ipFaceCoords[3](1,1) = a1;
ipFaceCoords[1](2,1) = a2; ipFaceCoords[3](2,1) = 0;
//
ipFaceCoords[1](0,2) = 0; ipFaceCoords[3](0,2) = a1;
ipFaceCoords[1](1,2) = a2; ipFaceCoords[3](1,2) = a2;
ipFaceCoords[1](2,2) = a1; ipFaceCoords[3](2,2) = 0;
// 2D integration points for faces
ipFace2DCoords(0,0) = a2; ipFace2DCoords(0,2) = a2;
ipFace2DCoords(1,0) = a2; ipFace2DCoords(1,2) = a1;
ipFace2DCoords(2,0) = a1; ipFace2DCoords(2,2) = a2;
//
ipFace2DCoords(0,1) = a1;
ipFace2DCoords(1,1) = a2;
ipFace2DCoords(2,1) = a2;
// Integration point weights
ipWeights = (1.0/6.0)/numIPs;
// Integration point face weights
ipFaceWeights = (1.0/2.0)/numFaceIPs;
break;
}
case GAUSS3: {
// Set number of IPs and face IPs
numIPs = 5;
numFaceIPs = 4;
// Size matrices and vectors thereof accordingly
ipCoords.resize(4,numIPs);
ipFaceCoords.assign(4,DENS_MAT(3,numFaceIPs));
ipFace2DCoords.resize(3,numFaceIPs);
ipWeights.reset(numIPs);
ipFaceWeights.reset(numFaceIPs);
double v1, v2, v3, a1, a2, a3;
/* These weights for calculating Gaussian Quadrature
* points are taken from the paper "Integration Points
* for Triangles and Tetrahedra Obtained from the
* Gaussian Quadrature Points for a Line" by K. S. Sunder
* and R. A. Cookson, Computers and Structures, Vol 21,
* No. 5, 1985. */
v1 = 1.0/4.0;
v2 = 1.0/2.0;
v3 = 1.0/6.0;
a1 = 1.0/3.0;
a2 = 3.0/5.0;
a3 = 1.0/5.0;
// Integration point coordinates
ipCoords.resize(4, numIPs);
ipCoords(0,0) = v1;
ipCoords(1,0) = v1;
ipCoords(2,0) = v1;
ipCoords(3,0) = v1;
//
ipCoords(0,1) = v3; ipCoords(0,3) = v3;
ipCoords(1,1) = v3; ipCoords(1,3) = v2;
ipCoords(2,1) = v3; ipCoords(2,3) = v3;
ipCoords(3,1) = v2; ipCoords(3,3) = v3;
//
ipCoords(0,2) = v2; ipCoords(0,4) = v3;
ipCoords(1,2) = v3; ipCoords(1,4) = v3;
ipCoords(2,2) = v3; ipCoords(2,4) = v2;
ipCoords(3,2) = v3; ipCoords(3,4) = v3;
// Integration points by face
// ...face 0 ...face 2
ipFaceCoords[0](0,0) = a1; ipFaceCoords[2](0,0) = a1;
ipFaceCoords[0](1,0) = a1; ipFaceCoords[2](1,0) = 0;
ipFaceCoords[0](2,0) = a1; ipFaceCoords[2](2,0) = a1;
//
ipFaceCoords[0](0,1) = a2; ipFaceCoords[2](0,1) = a3;
ipFaceCoords[0](1,1) = a3; ipFaceCoords[2](1,1) = 0;
ipFaceCoords[0](2,1) = a3; ipFaceCoords[2](2,1) = a3;
//
ipFaceCoords[0](0,2) = a3; ipFaceCoords[2](0,2) = a2;
ipFaceCoords[0](1,2) = a2; ipFaceCoords[2](1,2) = 0;
ipFaceCoords[0](2,2) = a3; ipFaceCoords[2](2,2) = a3;
//
ipFaceCoords[0](0,3) = a3; ipFaceCoords[2](0,3) = a3;
ipFaceCoords[0](1,3) = a3; ipFaceCoords[2](1,3) = 0;
ipFaceCoords[0](2,3) = a2; ipFaceCoords[2](2,3) = a2;
// ...face 1 ...face 3
ipFaceCoords[1](0,0) = 0; ipFaceCoords[3](0,0) = a1;
ipFaceCoords[1](1,0) = a1; ipFaceCoords[3](1,0) = a1;
ipFaceCoords[1](2,0) = a1; ipFaceCoords[3](2,0) = 0;
//
ipFaceCoords[1](0,1) = 0; ipFaceCoords[3](0,1) = a3;
ipFaceCoords[1](1,1) = a2; ipFaceCoords[3](1,1) = a3;
ipFaceCoords[1](2,1) = a3; ipFaceCoords[3](2,1) = 0;
//
ipFaceCoords[1](0,2) = 0; ipFaceCoords[3](0,2) = a3;
ipFaceCoords[1](1,2) = a3; ipFaceCoords[3](1,2) = a2;
ipFaceCoords[1](2,2) = a3; ipFaceCoords[3](2,2) = 0;
//
ipFaceCoords[1](0,3) = 0; ipFaceCoords[3](0,3) = a2;
ipFaceCoords[1](1,3) = a3; ipFaceCoords[3](1,3) = a3;
ipFaceCoords[1](2,3) = a2; ipFaceCoords[3](2,3) = 0;
// 2D integration points for faces
//
ipFace2DCoords(0,0) = a1; ipFace2DCoords(0,2) = a2;
ipFace2DCoords(1,0) = a1; ipFace2DCoords(1,2) = a3;
ipFace2DCoords(2,0) = a1; ipFace2DCoords(2,2) = a3;
//
ipFace2DCoords(0,1) = a3; ipFace2DCoords(0,3) = a3;
ipFace2DCoords(1,1) = a3; ipFace2DCoords(1,3) = a2;
ipFace2DCoords(2,1) = a2; ipFace2DCoords(2,3) = a3;
for (int i=0; i<numIPs; ++i) {
if (i < 1) ipWeights[i] = -4.0/5.0;
else ipWeights[i] = 9.0/20.0;
}
// Face integration point weights
for (int i=0; i<numFaceIPs; ++i) {
if (i < 1) ipFaceWeights[i] = -9.0/16.0;
else ipFaceWeights[i] = 25.0/48.0;
}
break;
}
// Error
default: {
throw ATC_Error("Unrecognized quadrature type "
"for element type TETRA.");
}
}
}
}
};
}; // namespace ATC
#endif // FE_QUADRATURE_H
c4science · Help