From 25514ecbe472cab87cacf2a92d7432c6dccac7a4 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 13 Nov 2004 23:57:28 +0000
Subject: [PATCH] - 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
---
Plugin/Integrate.cpp | 197 +++++++++++++++
Plugin/Integrate.h | 42 ++++
Plugin/ShapeFunctions.h | 524 ++++++++++++++++++++++++++++++++++++++++
3 files changed, 763 insertions(+)
create mode 100644 Plugin/Integrate.cpp
create mode 100644 Plugin/Integrate.h
create mode 100644 Plugin/ShapeFunctions.h
diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp
new file mode 100644
index 0000000000..f9f8112dd4
--- /dev/null
+++ b/Plugin/Integrate.cpp
@@ -0,0 +1,197 @@
+// $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;
+}
diff --git a/Plugin/Integrate.h b/Plugin/Integrate.h
new file mode 100644
index 0000000000..09fa2cb18b
--- /dev/null
+++ b/Plugin/Integrate.h
@@ -0,0 +1,42 @@
+#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
diff --git a/Plugin/ShapeFunctions.h b/Plugin/ShapeFunctions.h
new file mode 100644
index 0000000000..c4ebe19e8e
--- /dev/null
+++ b/Plugin/ShapeFunctions.h
@@ -0,0 +1,524 @@
+#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
--
GitLab