Skip to content
Snippets Groups Projects
Commit 25514ecb authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

- First draft (pretty much untested!) of new "Integrate" plugin to
  * integrate scalar fields over all the elements in a view
  * integrate the circulation of vector fields along line elements
  * integrate the flux of vector fields across surface elements

  Used with Plugin(DisplacementRaise) and Plugin(Evaluate) this
  permits for example to compute the area/volume of deformed
  configurations; and, with Plugin(CutPlane)+Plugin(Skin), the
  perimeter of deformed sections. Another interesting application is
  to use it on a vector field with Plugin(CutPlane), in order to
  compute fluxes across arbitrary cross-sections.

- Added "connectPoints" option to Plugin(CutParametric) so
  that we can feed its output to Plugin(Integrate)

- Added Normals and Tangents options to visualize the orientation of
  elements in post-processing views

- Added "swapOrientation" in Plugin(Transform) to change the
  orientation of the elements (in place) (+ moved the transformation
  routines from the view class into the plugin)

- fixed #defines in some of the plugin header files
parent 024e9517
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Plugin/Integrate.cpp 0 → 100644
+ 197
− 0
View file @ 25514ecb
// $Id: Integrate.cpp,v 1.1 2004-11-13 23:57:28 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
// 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.
//
// This program 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 this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "Plugin.h"
#include "Integrate.h"
#include "List.h"
#include "Views.h"
#include "Context.h"
#include "Numeric.h"
#include "ShapeFunctions.h"
extern Context_T CTX;
StringXNumber IntegrateOptions_Number[] = {
{GMSH_FULLRC, "iView", NULL, -1.}
};
extern "C"
{
GMSH_Plugin *GMSH_RegisterIntegratePlugin()
{
return new GMSH_IntegratePlugin();
}
}
GMSH_IntegratePlugin::GMSH_IntegratePlugin()
{
;
}
void GMSH_IntegratePlugin::getName(char *name) const
{
strcpy(name, "Integrate");
}
void GMSH_IntegratePlugin::getInfos(char *author, char *copyright,
char *help_text) const
{
strcpy(author, "C. Geuzaine (geuz@geuz.org)");
strcpy(copyright, "DGR (www.multiphysics.com)");
strcpy(help_text,
"Plugin(Integrate) integrates scalar fields over\n"
"all the elements in the view `iView', as well as\n"
"the circulation/flux of vector fields over\n"
"line/surface elements. If `iView' < 0, the plugin\n"
"is run on the current view.\n"
"\n"
"Plugin(Integrate) creates one new view.\n");
}
int GMSH_IntegratePlugin::getNbOptions() const
{
return sizeof(IntegrateOptions_Number) / sizeof(StringXNumber);
}
StringXNumber *GMSH_IntegratePlugin::getOption(int iopt)
{
return &IntegrateOptions_Number[iopt];
}
void GMSH_IntegratePlugin::catchErrorMessage(char *errorMessage) const
{
strcpy(errorMessage, "Integrate failed...");
}
static double integrate(int nbList, List_T *list, int dim,
int nbNod, int nbComp, int step)
{
if(!nbList) return 0.;
double res = 0.;
int nb = List_Nbr(list) / nbList;
for(int i = 0; i < List_Nbr(list); i += nb) {
double *x = (double *)List_Pointer_Fast(list, i);
double *y = (double *)List_Pointer_Fast(list, i + nbNod);
double *z = (double *)List_Pointer_Fast(list, i + 2 * nbNod);
double *v = (double *)List_Pointer_Fast(list, i + 3 * nbNod +
nbNod * nbComp * step);
if(dim == 0){
if(nbNod == 1){
point p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
}
}
else if(dim == 1){
if(nbNod == 2){
line l(x, y, z);
if(nbComp == 1) res += l.integrate(v);
else if(nbComp == 3) res += l.integrateCirculation(v);
}
}
else if(dim == 2){
if(nbNod == 3){
triangle t(x, y, z);
if(nbComp == 1) res += t.integrate(v);
else if(nbComp == 3) res += t.integrateFlux(v);
}
else if(nbNod == 4){
quadrangle q(x, y, z);
if(nbComp == 1) res += q.integrate(v);
else if(nbComp == 3) res += q.integrateFlux(v);
}
}
else if(dim == 3){
if(nbNod == 4){
tetrahedron t(x, y, z);
if(nbComp == 1) res += t.integrate(v);
}
else if(nbNod == 8){
hexahedron h(x, y, z);
if(nbComp == 1) res += h.integrate(v);
}
else if(nbNod == 6){
prism p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
}
else if(nbNod == 5){
pyramid p(x, y, z);
if(nbComp == 1) res += p.integrate(v);
}
}
}
return res;
}
Post_View *GMSH_IntegratePlugin::execute(Post_View * v)
{
int iView = (int)IntegrateOptions_Number[0].def;
if(iView < 0)
iView = v ? v->Index : 0;
if(!List_Pointer_Test(CTX.post.list, iView)) {
Msg(GERROR, "View[%d] does not exist", iView);
return v;
}
Post_View *v2 = BeginView(1);
Post_View *v1 = (Post_View*)List_Pointer(CTX.post.list, iView);
double x = (v1->BBox[0]+v1->BBox[1])/2.;
double y = (v1->BBox[2]+v1->BBox[3])/2.;
double z = (v1->BBox[4]+v1->BBox[5])/2.;
List_Add(v2->SP, &x);
List_Add(v2->SP, &y);
List_Add(v2->SP, &z);
for(int ts = 0; ts < v1->NbTimeStep; ts++){
double val = 0;
// scalar fields
val += integrate(v1->NbSP, v1->SP, 0, 1, 1, ts);
val += integrate(v1->NbSL, v1->SL, 1, 2, 1, ts);
val += integrate(v1->NbST, v1->ST, 2, 3, 1, ts);
val += integrate(v1->NbSQ, v1->SQ, 2, 4, 1, ts);
val += integrate(v1->NbSS, v1->SS, 3, 4, 1, ts);
val += integrate(v1->NbSH, v1->SH, 3, 8, 1, ts);
val += integrate(v1->NbSI, v1->SI, 3, 6, 1, ts);
val += integrate(v1->NbSY, v1->SY, 3, 5, 1, ts);
// circulations
val += integrate(v1->NbVL, v1->VL, 1, 2, 3, ts);
// fluxes
val += integrate(v1->NbVT, v1->VT, 2, 3, 3, ts);
val += integrate(v1->NbVQ, v1->VQ, 2, 4, 3, ts);
Msg(INFO, "Step %d: integral = %.16g", ts, val);
List_Add(v2->SP, &val);
}
v2->NbSP = 1;
v2->IntervalsType = DRAW_POST_NUMERIC;
// copy time data
for(int i = 0; i < List_Nbr(v1->Time); i++)
List_Add(v2->Time, List_Pointer(v1->Time, i));
// finalize
char name[1024], filename[1024];
sprintf(name, "%s_Integrate", v1->Name);
sprintf(filename, "%s_Integrate.pos", v1->Name);
EndView(v2, 1, filename, name);
return v2;
}
Plugin/Integrate.h 0 → 100644
+ 42
− 0
View file @ 25514ecb
#ifndef _INTEGRATE_H_
#define _INTEGRATE_H_
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
// 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.
//
// This program 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 this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "Plugin.h"
extern "C"
{
GMSH_Plugin *GMSH_RegisterIntegratePlugin();
}
class GMSH_IntegratePlugin : public GMSH_Post_Plugin
{
public:
GMSH_IntegratePlugin();
void getName(char *name) const;
void getInfos(char *author, char *copyright, char *help_text) const;
void catchErrorMessage(char *errorMessage) const;
int getNbOptions() const;
StringXNumber* getOption(int iopt);
Post_View *execute(Post_View *);
};
#endif
Plugin/ShapeFunctions.h 0 → 100644
+ 524
− 0
View file @ 25514ecb
#ifndef _SHAPE_FUNCTIONS_H_
#define _SHAPE_FUNCTIONS_H_
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
// 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.
//
// This program 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 this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "Numeric.h"
class element{
protected:
double *_x, *_y, *_z;
public:
element(double *x, double *y, double *z) : _x(x), _y(y), _z(z) {}
virtual ~element(){}
virtual int getDimension() = 0;
virtual int getNumNodes() = 0;
virtual int getNumGaussPoints() = 0;
virtual void getGaussPoint(int num, double &u, double &v, double &w, double &weight) = 0;
virtual void getShapeFunction(int num, double u, double v, double w, double &s) = 0;
virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]) = 0;
double getJacobian(double u, double v, double w, double jac[3][3])
{
jac[0][0] = jac[0][1] = jac[0][2] = 0.;
jac[1][0] = jac[1][1] = jac[1][2] = 0.;
jac[2][0] = jac[2][1] = jac[2][2] = 0.;
double s[3];
switch(getDimension()){
case 3 :
for(int i = 0; i < getNumNodes(); i++) {
getGradShapeFunction(i, u, v, w, s);
jac[0][0] += _x[i] * s[0]; jac[0][1] += _y[i] * s[0]; jac[0][2] += _z[i] * s[0];
jac[1][0] += _x[i] * s[1]; jac[1][1] += _y[i] * s[1]; jac[1][2] += _z[i] * s[1];
jac[2][0] += _x[i] * s[2]; jac[2][1] += _y[i] * s[2]; jac[2][2] += _z[i] * s[2];
}
return
jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2];
case 2 :
for(int i = 0; i < getNumNodes(); i++) {
getGradShapeFunction(i, u, v, w, s);
jac[0][0] += _x[i] * s[0]; jac[0][1] += _y[i] * s[0]; jac[0][2] += _z[i] * s[0];
jac[1][0] += _x[i] * s[1]; jac[1][1] += _y[i] * s[1]; jac[1][2] += _z[i] * s[1];
}
return sqrt(DSQR(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
DSQR(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
DSQR(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
case 1:
for(int i = 0; i < getNumNodes(); i++) {
getGradShapeFunction(i, u, v, w, s);
jac[0][0] += _x[i] * s[0]; jac[0][1] += _y[i] * s[0]; jac[0][2] += _z[i] * s[0];
}
return sqrt(DSQR(jac[0][0])+DSQR(jac[0][1])+DSQR(jac[0][2]));
default:
return 1.;
}
}
double interpolate(double val[], double u, double v, double w)
{
double sum = 0;
for(int i = 0; i < getNumNodes(); i++){
double s;
getShapeFunction(i, u, v, w, s);
sum += val[i] * s;
}
return sum;
}
double integrate(double val[])
{
double sum = 0;
for(int i = 0; i < getNumGaussPoints(); i++){
double u, v, w, weight, jac[3][3];
getGaussPoint(i, u, v, w, weight);
double det = getJacobian(u, v, w, jac);
double d = interpolate(val, u, v, w);
sum += d * weight * det;
}
return sum;
}
};
class point : public element{
public:
point(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 0; }
inline int getNumNodes(){ return 1; }
inline int getNumGaussPoints(){ return 1; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
u = v = w = 0.;
weight = 1.;
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = 1.; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
s[0] = s[1] = s[2] = 0.;
}
};
class line : public element{
public:
line(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 1; }
inline int getNumNodes(){ return 2; }
inline int getNumGaussPoints(){ return 1; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
if(num < 0 || num > 0) return;
u = v = w = 0.;
weight = 2.;
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = 0.5 * (1.-u); break;
case 1 : s = 0.5 * (1.+u); break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -0.5; s[1] = 0.; s[2] = 0.; break;
case 1 : s[0] = 0.5; s[1] = 0.; s[2] = 0.; break;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
double integrateCirculation(double val[])
{
double t[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
norme(t);
double v[3];
for(int i = 0; i < 3; i++){
double tmp[2];
tmp[0] = val[i];
tmp[1] = val[3+i];
v[i] = integrate(tmp);
}
double d;
prosca(t, v, &d);
return d;
}
};
class triangle : public element{
public:
triangle(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 2; }
inline int getNumNodes(){ return 3; }
inline int getNumGaussPoints(){ return 3; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u3[3] = {0.16666666666666,0.66666666666666,0.16666666666666};
static double v3[3] = {0.16666666666666,0.16666666666666,0.66666666666666};
static double p3[3] = {0.16666666666666,0.16666666666666,0.16666666666666};
if(num < 0 || num > 2) return;
u = u3[num];
v = v3[num];
w = 0.;
weight = p3[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num){
case 0 : s = 1.-u-v; break;
case 1 : s = u ; break;
case 2 : s = v; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -1.; s[1] = -1.; s[2] = 0.; break;
case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
double integrateFlux(double val[])
{
double t1[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
double t2[3] = {_x[2]-_x[0], _y[2]-_y[0], _z[2]-_z[0]};
double n[3];
prodve(t1, t2, n);
norme(n);
double v[3];
for(int i = 0; i < 3; i++){
double tmp[3];
tmp[0] = val[i];
tmp[1] = val[3+i];
tmp[2] = val[6+i];
v[i] = integrate(tmp);
}
double d;
prosca(n, v, &d);
return d;
}
};
class quadrangle : public element{
public:
quadrangle(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 2; }
inline int getNumNodes(){ return 4; }
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u4[4] = {0.577350269189,-0.577350269189,0.577350269189,-0.577350269189};
static double v4[4] = {0.577350269189,0.577350269189,-0.577350269189,-0.577350269189};
static double p4[4] = {1.,1.,1.,1.};
if(num < 0 || num > 3) return;
u = u4[num];
v = v4[num];
w = 0.;
weight = p4[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = 0.25 * (1.-u) * (1.-v); break;
case 1 : s = 0.25 * (1.+u) * (1.-v); break;
case 2 : s = 0.25 * (1.+u) * (1.+v); break;
case 3 : s = 0.25 * (1.-u) * (1.+v); break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -0.25 * (1.-v); s[1] = -0.25 * (1.-u); s[2] = 0.; break;
case 1 : s[0] = 0.25 * (1.-v); s[1] = -0.25 * (1.+u); s[2] = 0.; break;
case 2 : s[0] = 0.25 * (1.+v); s[1] = 0.25 * (1.+u); s[2] = 0.; break;
case 3 : s[0] = -0.25 * (1.+v); s[1] = 0.25 * (1.-u); s[2] = 0.; break;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
double integrateFlux(double val[])
{
double t1[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
double t2[3] = {_x[2]-_x[0], _y[2]-_y[0], _z[2]-_z[0]};
double n[3];
prodve(t1, t2, n);
norme(n);
double v[3];
for(int i = 0; i < 3; i++){
double tmp[4];
tmp[0] = val[i];
tmp[1] = val[3+i];
tmp[2] = val[6+i];
tmp[3] = val[9+i];
v[i] = integrate(tmp);
}
double d;
prosca(n, v, &d);
return d;
}
};
class tetrahedron : public element{
public:
tetrahedron(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 4; }
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u4[4] = {0.138196601125,0.138196601125,0.138196601125,0.585410196625};
static double v4[4] = {0.138196601125,0.138196601125,0.585410196625,0.138196601125};
static double w4[4] = {0.138196601125,0.585410196625,0.138196601125,0.138196601125};
static double p4[4] = {0.0416666666667,0.0416666666667,0.0416666666667,0.0416666666667};
if(num < 0 || num > 3) return;
u = u4[num];
v = v4[num];
w = w4[num];
weight = p4[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = 1.-u-v-w; break;
case 1 : s = u ; break;
case 2 : s = v ; break;
case 3 : s = w; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -1.; s[1] = -1.; s[2] = -1.; break;
case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
case 3 : s[0] = 0.; s[1] = 0.; s[2] = 1.; break;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
};
class hexahedron : public element{
public:
hexahedron(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 8; }
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u6[6] = { 0.40824826, 0.40824826, -0.40824826,
-0.40824826, -0.81649658, 0.81649658};
static double v6[6] = { 0.70710678, -0.70710678, 0.70710678,
-0.70710678, 0., 0.};
static double w6[6] = {-0.57735027, -0.57735027, 0.57735027,
0.57735027, -0.57735027, 0.57735027};
static double p6[6] = { 1.3333333333, 1.3333333333, 1.3333333333,
1.3333333333, 1.3333333333, 1.3333333333};
if(num < 0 || num > 5) return;
u = u6[num];
v = v6[num];
w = w6[num];
weight = p6[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = (1.-u) * (1.-v) * (1.-w) * 0.125; break;
case 1 : s = (1.+u) * (1.-v) * (1.-w) * 0.125; break;
case 2 : s = (1.+u) * (1.+v) * (1.-w) * 0.125; break;
case 3 : s = (1.-u) * (1.+v) * (1.-w) * 0.125; break;
case 4 : s = (1.-u) * (1.-v) * (1.+w) * 0.125; break;
case 5 : s = (1.+u) * (1.-v) * (1.+w) * 0.125; break;
case 6 : s = (1.+u) * (1.+v) * (1.+w) * 0.125; break;
case 7 : s = (1.-u) * (1.+v) * (1.+w) * 0.125; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -0.125 * (1.-v) * (1.-w);
s[1] = -0.125 * (1.-u) * (1.-w);
s[2] = -0.125 * (1.-u) * (1.-v); break;
case 1 : s[0] = 0.125 * (1.-v) * (1.-w);
s[1] = -0.125 * (1.+u) * (1.-w);
s[2] = -0.125 * (1.+u) * (1.-v); break;
case 2 : s[0] = 0.125 * (1.+v) * (1.-w);
s[1] = 0.125 * (1.+u) * (1.-w);
s[2] = -0.125 * (1.+u) * (1.+v); break;
case 3 : s[0] = -0.125 * (1.+v) * (1.-w);
s[1] = 0.125 * (1.-u) * (1.-w);
s[2] = -0.125 * (1.-u) * (1.+v); break;
case 4 : s[0] = -0.125 * (1.-v) * (1.+w);
s[1] = -0.125 * (1.-u) * (1.+w);
s[2] = 0.125 * (1.-u) * (1.-v); break;
case 5 : s[0] = 0.125 * (1.-v) * (1.+w);
s[1] = -0.125 * (1.+u) * (1.+w);
s[2] = 0.125 * (1.+u) * (1.-v); break;
case 6 : s[0] = 0.125 * (1.+v) * (1.+w);
s[1] = 0.125 * (1.+u) * (1.+w);
s[2] = 0.125 * (1.+u) * (1.+v); break;
case 7 : s[0] = -0.125 * (1.+v) * (1.+w);
s[1] = 0.125 * (1.-u) * (1.+w);
s[2] = 0.125 * (1.-u) * (1.+v); break;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
};
class prism : public element{
public:
prism(double *x, double *y, double *z) : element(x, y, z) {};
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 6; }
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u6[6] = {0.166666666666666, 0.333333333333333, 0.166666666666666,
0.166666666666666, 0.333333333333333, 0.166666666666666};
static double v6[6] = {0.166666666666666, 0.166666666666666, 0.333333333333333,
0.166666666666666, 0.166666666666666, 0.333333333333333};
static double w6[6] = {-0.577350269189, -0.577350269189, -0.577350269189,
0.577350269189, 0.577350269189, 0.577350269189};
static double p6[6] = {0.166666666666666,0.166666666666666,0.166666666666666,
0.166666666666666,0.166666666666666,0.166666666666666,};
if(num < 0 || num > 5) return;
u = u6[num];
v = v6[num];
w = w6[num];
weight = p6[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
switch(num) {
case 0 : s = (1.-u-v) * (1.-w) * 0.5; break;
case 1 : s = u * (1.-w) * 0.5; break;
case 2 : s = v * (1.-w) * 0.5; break;
case 3 : s = (1.-u-v) * (1.+w) * 0.5; break;
case 4 : s = u * (1.+w) * 0.5; break;
case 5 : s = v * (1.+w) * 0.5; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
case 0 : s[0] = -0.5 * (1.-w) ;
s[1] = -0.5 * (1.-w) ;
s[2] = -0.5 * (1.-u-v) ; break ;
case 1 : s[0] = 0.5 * (1.-w) ;
s[1] = 0. ;
s[2] = -0.5 * u ; break ;
case 2 : s[0] = 0. ;
s[1] = 0.5 * (1.-w) ;
s[2] = -0.5 * v ; break ;
case 3 : s[0] = -0.5 * (1.+w) ;
s[1] = -0.5 * (1.+w) ;
s[2] = 0.5 * (1.-u-v) ; break ;
case 4 : s[0] = 0.5 * (1.+w) ;
s[1] = 0. ;
s[2] = 0.5 * u ; break ;
case 5 : s[0] = 0. ;
s[1] = 0.5 * (1.+w) ;
s[2] = 0.5 * v ; break ;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
};
class pyramid : public element{
public:
pyramid(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 5; }
inline int getNumGaussPoints(){ return 8; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
static double u8[8] = {0.3595161057791018,0.09633205020967324,
0.3595161057791018,0.09633205020967324,
0.6920507403468987,0.1854344369976602,
0.6920507403468987,0.1854344369976602};
static double v8[8] = {0.3595161057791018,0.3595161057791018,
0.09633205020967324,0.09633205020967324,
0.6920507403468987,0.6920507403468987,
0.1854344369976602,0.1854344369976602};
static double w8[8] = {0.544151844011225,0.544151844011225,
0.544151844011225,0.544151844011225,
0.122514822655441,0.122514822655441,
0.122514822655441,0.122514822655441};
static double p8[8] = {0.02519647051995625,0.02519647051995625,
0.02519647051995625,0.02519647051995625,
0.058136862813377,0.058136862813377,
0.058136862813377,0.058136862813377};
if(num < 0 || num > 7) return;
u = u8[num];
v = v8[num];
w = w8[num];
weight = p8[num];
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
double r;
if(w != 1. && num != 4) r = u*v*w / (1.-w);
else r = 0.;
switch(num) {
case 0 : s = 0.25 * ((1.-u) * (1.-v) - w + r); break;
case 1 : s = 0.25 * ((1.+u) * (1.-v) - w - r); break;
case 2 : s = 0.25 * ((1.+u) * (1.+v) - w + r); break;
case 3 : s = 0.25 * ((1.-u) * (1.+v) - w - r); break;
case 4 : s = w; break;
default : s = 0.; break;
}
}
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
if(w == 1. && num != 4) {
s[0] = 0.25;
s[1] = 0.25;
s[2] = -0.25;
}
else{
switch(num) {
case 0 : s[0] = 0.25 * ( -(1.-v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( -(1.-u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/DSQR(1.-w) ) ; break ;
case 1 : s[0] = 0.25 * ( (1.-v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( -(1.+u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/DSQR(1.-w) ) ; break ;
case 2 : s[0] = 0.25 * ( (1.+v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( (1.+u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/DSQR(1.-w) ) ; break ;
case 3 : s[0] = 0.25 * ( -(1.+v) + v*w/(1.-w) ) ;
s[1] = 0.25 * ( (1.-u) + u*w/(1.-w) ) ;
s[2] = 0.25 * ( -1. + u*v/DSQR(1.-w) ) ; break ;
case 4 : s[0] = 0. ;
s[1] = 0. ;
s[2] = 1. ; break ;
default : s[0] = s[1] = s[2] = 0.; break;
}
}
}
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment