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

add an elementFactory to simplify the tests + new "getNode" fcts
parent eda9d43a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Plugin/Integrate.cpp
+ 17
− 66
View file @ f8ccc36a
// $Id: Integrate.cpp,v 1.10 2004-11-26 16:46:52 geuzaine Exp $
// $Id: Integrate.cpp,v 1.11 2004-12-13 04:00:18 geuzaine Exp $
//
// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
//
......@@ -99,75 +99,26 @@ static double integrate(int nbList, List_T *list, int dim,
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){
if(!levelsetPositive) res += p.integrate(v);
else res += p.integrateLevelsetPositive(v);
}
}
elementFactory factory;
element *element = factory.create(nbNod, dim, x, y, z);
if(!element){
Msg(GERROR, "Unknown type of element (dim=%d, nbNod=%d)", dim, nbNod);
return 0.;
}
else if(dim == 1){
if(nbNod == 2){
line l(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += l.integrate(v);
else res += l.integrateLevelsetPositive(v);
}
else if(nbComp == 3) res += l.integrateCirculation(v);
}
if(nbComp == 1){
if(!levelsetPositive)
res += element->integrate(v);
else
res += element->integrateLevelsetPositive(v);
}
else if(dim == 2){
if(nbNod == 3){
triangle t(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += t.integrate(v);
else res += t.integrateLevelsetPositive(v);
}
else if(nbComp == 3) res += t.integrateFlux(v);
}
else if(nbNod == 4){
quadrangle q(x, y, z);
if(nbComp == 1) {
if(!levelsetPositive) res += q.integrate(v);
else res += q.integrateLevelsetPositive(v);
}
else if(nbComp == 3) res += q.integrateFlux(v);
}
}
else if(dim == 3){
if(nbNod == 4){
tetrahedron t(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += t.integrate(v);
else res += t.integrateLevelsetPositive(v);
}
}
else if(nbNod == 8){
hexahedron h(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += h.integrate(v);
else res += h.integrateLevelsetPositive(v);
}
}
else if(nbNod == 6){
prism p(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += p.integrate(v);
else res += p.integrateLevelsetPositive(v);
}
}
else if(nbNod == 5){
pyramid p(x, y, z);
if(nbComp == 1){
if(!levelsetPositive) res += p.integrate(v);
else res += p.integrateLevelsetPositive(v);
}
}
else if(nbComp == 3){
if(dim == 1)
res += element->integrateCirculation(v);
else if(dim == 2)
res += element->integrateFlux(v);
}
delete element;
}
// printf("integration res = %22.15E\n",res);
return res;
}
......
Plugin/ShapeFunctions.h
+ 142
− 3
View file @ f8ccc36a
......@@ -30,6 +30,7 @@ public:
virtual ~element(){}
virtual int getDimension() = 0;
virtual int getNumNodes() = 0;
virtual void getNode(int num, double &u, double &v, double &w) = 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;
......@@ -71,24 +72,52 @@ public:
return 1.;
}
}
double interpolate(double val[], double u, double v, double w)
double getInverseJacobian(double jac[3][3], double invjac[3][3])
{
// maybe we should use the classic approach to define always
// non-singular jacobians, so that we can simply invert them with
// inv3x3
Msg(GERROR, "getInverseJacibian not done yet");
}
double interpolate(double val[], double u, double v, double w, int stride=1)
{
double sum = 0;
int j = 0;
for(int i = 0; i < getNumNodes(); i++){
double s;
getShapeFunction(i, u, v, w, s);
sum += val[i] * s;
j += stride;
}
return sum;
}
double integrate(double val[])
void interpolateGrad(double val[], double u, double v, double w, double f[3], int stride=1)
{
double fu[3] = {0., 0., 0.};
int j = 0;
for(int i = 0; i < getNumNodes(); i++){
double s[3];
getGradShapeFunction(i, u, v, w, s);
fu[0] += val[j] * s[0];
fu[1] += val[j] * s[1];
fu[2] += val[j] * s[2];
j += stride;
}
double jac[3][3], invjac[3][3];
getJacobian(u, v, w, jac);
getInverseJacobian(jac, invjac);
f[0] = fu[0] * invjac[0][0] + fu[1] * invjac[0][1] + fu[2] * invjac[0][2];
f[1] = fu[0] * invjac[1][0] + fu[1] * invjac[1][1] + fu[2] * invjac[1][2];
f[2] = fu[0] * invjac[2][0] + fu[1] * invjac[2][1] + fu[2] * invjac[2][2];
}
double integrate(double val[], int stride=1)
{
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);
double d = interpolate(val, u, v, w, stride);
sum += d * weight * det;
}
return sum;
......@@ -107,6 +136,16 @@ public:
res = area * (1 - sum/sumabs) * 0.5 ;
return res;
}
virtual double integrateCirculation(double val[])
{
Msg(GERROR, "integrateCirculation not available for this element");
return 0.;
}
virtual double integrateFlux(double val[])
{
Msg(GERROR, "integrateFlux not available for this element");
return 0.;
}
};
class point : public element{
......@@ -114,6 +153,10 @@ public:
point(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 0; }
inline int getNumNodes(){ return 1; }
void getNode(int num, double &u, double &v, double &w)
{
u = v = w = 0.;
}
inline int getNumGaussPoints(){ return 1; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -138,6 +181,15 @@ public:
line(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 1; }
inline int getNumNodes(){ return 2; }
void getNode(int num, double &u, double &v, double &w)
{
v = w = 0.;
switch(num) {
case 0 : u = -1.; break;
case 1 : u = 1.; break;
default: u = 0.; break;
}
}
inline int getNumGaussPoints(){ return 1; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -183,6 +235,16 @@ public:
triangle(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 2; }
inline int getNumNodes(){ return 3; }
void getNode(int num, double &u, double &v, double &w)
{
w = 0.;
switch(num) {
case 0 : u = 0.; v = 0.; break;
case 1 : u = 1.; v = 0.; break;
case 2 : u = 0.; v = 1.; break;
default: u = 0.; v = 0.; break;
}
}
inline int getNumGaussPoints(){ return 3; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -239,6 +301,17 @@ public:
quadrangle(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 2; }
inline int getNumNodes(){ return 4; }
void getNode(int num, double &u, double &v, double &w)
{
w = 0.;
switch(num) {
case 0 : u = -1.; v = -1.; break;
case 1 : u = 1.; v = -1.; break;
case 2 : u = 1.; v = 1.; break;
case 3 : u = -1.; v = 1.; break;
default: u = 0.; v = 0.; break;
}
}
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -298,6 +371,16 @@ public:
tetrahedron(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 4; }
void getNode(int num, double &u, double &v, double &w)
{
switch(num) {
case 0 : u = 0.; v = 0.; w = 0.; break;
case 1 : u = 1.; v = 0.; w = 0.; break;
case 2 : u = 0.; v = 1.; w = 0.; break;
case 3 : u = 0.; v = 0.; w = 1.; break;
default: u = 0.; v = 0.; w = 0.; break;
}
}
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -338,6 +421,20 @@ public:
hexahedron(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 8; }
void getNode(int num, double &u, double &v, double &w)
{
switch(num) {
case 0 : u = -1.; v = -1.; w = -1.; break;
case 1 : u = 1.; v = -1.; w = -1.; break;
case 2 : u = 1.; v = 1.; w = -1.; break;
case 3 : u = -1.; v = 1.; w = -1.; break;
case 4 : u = -1.; v = -1.; w = 1.; break;
case 5 : u = 1.; v = -1.; w = 1.; break;
case 6 : u = 1.; v = 1.; w = 1.; break;
case 7 : u = -1.; v = 1.; w = 1.; break;
default: u = 0.; v = 0.; w = 0.; break;
}
}
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -406,6 +503,18 @@ public:
prism(double *x, double *y, double *z) : element(x, y, z) {};
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 6; }
void getNode(int num, double &u, double &v, double &w)
{
switch(num) {
case 0 : u = 0.; v = 0.; w = -1.; break;
case 1 : u = 1.; v = 0.; w = -1.; break;
case 2 : u = 0.; v = 1.; w = -1.; break;
case 3 : u = 0.; v = 0.; w = 1.; break;
case 4 : u = 1.; v = 0.; w = 1.; break;
case 5 : u = 0.; v = 1.; w = 1.; break;
default: u = 0.; v = 0.; w = 0.; break;
}
}
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -466,6 +575,17 @@ public:
pyramid(double *x, double *y, double *z) : element(x, y, z) {}
inline int getDimension(){ return 3; }
inline int getNumNodes(){ return 5; }
void getNode(int num, double &u, double &v, double &w)
{
switch(num) {
case 0 : u = -1.; v = -1.; w = 0.; break;
case 1 : u = 1.; v = -1.; w = 0.; break;
case 2 : u = 1.; v = 1.; w = 0.; break;
case 3 : u = -1.; v = 1.; w = 0.; break;
case 4 : u = 0.; v = 0.; w = 1.; break;
default: u = 0.; v = 0.; w = 0.; break;
}
}
inline int getNumGaussPoints(){ return 8; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
......@@ -535,4 +655,23 @@ public:
}
};
class elementFactory{
public:
element* create(int numNodes, int dimension, double *x, double *y, double *z){
switch(dimension){
case 0: return new point(x, y, z);
case 1: return new line(x, y, z);
case 2:
if(numNodes == 4) return new quadrangle(x, y, z);
else return new triangle(x, y, z);
case 3:
if(numNodes == 8) return new hexahedron(x, y, z);
else if(numNodes == 6) return new prism(x, y, z);
else if(numNodes == 5) return new pyramid(x, y, z);
else return new tetrahedron(x, y, z);
default: return NULL;
}
}
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment