From 59c8dd78bec98f5ec85c480e204660d504641356 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 24 Jun 2012 11:58:00 +0000
Subject: [PATCH] work on far field
---
Plugin/Integrate.cpp | 23 ++---
Plugin/NearToFarField.cpp | 209 ++++++++++++++++++++++++++------------
Plugin/NearToFarField.h | 15 ++-
Post/shapeFunctions.h | 101 +++++++++---------
4 files changed, 217 insertions(+), 131 deletions(-)
diff --git a/Plugin/Integrate.cpp b/Plugin/Integrate.cpp
index 5746451d7d..e4a2f1d222 100644
--- a/Plugin/Integrate.cpp
+++ b/Plugin/Integrate.cpp
@@ -27,7 +27,7 @@ std::string GMSH_IntegratePlugin::getHelp() const
"as the circulation/flux of vector fields over "
"line/surface elements.\n\n"
"If `View' < 0, the plugin is run on the current view.\n\n"
- "Plugin(Integrate) creates one new view.";
+ "Plugin(Integrate) creates one new view."
"If `OverTime' = 1 , the plugin integrates the scalar view over time instead of over space.\n\n"
"Plugin(Integrate) creates one new view.";
}
@@ -46,15 +46,15 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
{
int iView = (int)IntegrateOptions_Number[0].def;
int overTime = (int)IntegrateOptions_Number[1].def;
-
+
PView *v1 = getView(iView, v);
if(!v1) return v;
-
+
PViewData *data1 = getPossiblyAdaptiveData(v1);
PView *v2 = new PView();
PViewDataList *data2 = getDataList(v2);
-
- if (overTime == -1) {
+
+ if (overTime == -1) {
double x = data1->getBoundingBox().center().x();
double y = data1->getBoundingBox().center().y();
double z = data1->getBoundingBox().center().z();
@@ -83,7 +83,7 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
if(numNodes == 1){
simpleSum = true;
res += val[0];
- for(int comp = 0; comp < numComp; comp++)
+ for(int comp = 0; comp < numComp; comp++)
resv[comp] += val[comp];
}
else{
@@ -101,8 +101,8 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
}
}
if(simpleSum)
- Msg::Info("Step %d: sum = %g %g %g %g %g %g %g %g %g", step, resv[0],
- resv[1], resv[2], resv[3], resv[4], resv[5], resv[6], resv[7],
+ Msg::Info("Step %d: sum = %g %g %g %g %g %g %g %g %g", step, resv[0],
+ resv[1], resv[2], resv[3], resv[4], resv[5], resv[6], resv[7],
resv[8]);
else
Msg::Info("Step %d: integral = %.16g", step, res);
@@ -110,7 +110,7 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
}
data2->NbSP = 1;
v2->getOptions()->intervalsType = PViewOptions::Numeric;
-
+
for(int i = 0; i < data1->getNumTimeSteps(); i++){
double time = data1->getTime(i);
data2->Time.push_back(time);
@@ -136,12 +136,11 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
for(int nod = 0; nod < numNodes; nod++) out->push_back(x[nod]);
for(int nod = 0; nod < numNodes; nod++) out->push_back(y[nod]);
for(int nod = 0; nod < numNodes; nod++) out->push_back(z[nod]);
-
+
double time = 0.0;
std::vector<double> timeIntegral(numNodes, 0.);
for(int step = timeBeg; step < timeEnd; step++){
if(!data1->hasTimeStep(step)) continue;
- double val;
double newTime = data1->getTime(step);
double dt = newTime - time;
time = newTime;
@@ -160,6 +159,6 @@ PView *GMSH_IntegratePlugin::execute(PView * v)
data2->setName(data1->getName() + "_Integrate");
data2->setFileName(data1->getName() + "_Integrate.pos");
data2->finalize();
-
+
return v2;
}
diff --git a/Plugin/NearToFarField.cpp b/Plugin/NearToFarField.cpp
index 4f6a08375d..0c532e84cd 100644
--- a/Plugin/NearToFarField.cpp
+++ b/Plugin/NearToFarField.cpp
@@ -2,18 +2,27 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Ruth Sabariego
+//
+#include <complex>
#include "NearToFarField.h"
StringXNumber NearToFarFieldOptions_Number[] = {
- {GMSH_FULLRC, "Wavenumber", NULL, 1.},
- {GMSH_FULLRC, "FarDistance", NULL, 1.},
- {GMSH_FULLRC, "NumPointsPhi", NULL, 120},
- {GMSH_FULLRC, "NumPointsTheta", NULL, 60},
- {GMSH_FULLRC, "EView", NULL, 0},
- {GMSH_FULLRC, "HView", NULL, 1},
- {GMSH_FULLRC, "Normalize", NULL, 1},
- {GMSH_FULLRC, "dB", NULL, 1},
+ {GMSH_FULLRC, "Wavenumber", NULL, 1.},
+ {GMSH_FULLRC, "PhiStart", NULL, 0.},
+ {GMSH_FULLRC, "PhiEnd", NULL, 2. * M_PI},
+ {GMSH_FULLRC, "NumPointsPhi", NULL, 60},
+ {GMSH_FULLRC, "ThetaStart", NULL, 0.},
+ {GMSH_FULLRC, "ThetaEnd", NULL, M_PI},
+ {GMSH_FULLRC, "NumPointsTheta", NULL, 30},
+ {GMSH_FULLRC, "EView", NULL, 0},
+ {GMSH_FULLRC, "HView", NULL, 1},
+ {GMSH_FULLRC, "Normalize", NULL, 1},
+ {GMSH_FULLRC, "dB", NULL, 1},
+ {GMSH_FULLRC, "NegativeTime", NULL, 0.},
};
extern "C"
@@ -27,12 +36,16 @@ extern "C"
std::string GMSH_NearToFarFieldPlugin::getHelp() const
{
return "Plugin(NearToFarField) computes the far field pattern "
- "from the near electric and magnetic fields on a surface "
+ "from the near electric E and magnetic H fields on a surface "
"enclosing the radiating device (antenna).\n\n"
- "Parameters: the wavenumber, the far field distance (radius) and "
- "angular discretisation, i.e. the number of divisions for "
- "phi in [0, 2*Pi] and theta in [0, Pi].\n\n"
- "If `View' < 0, the plugin is run on the current view.\n\n"
+ "Parameters: the wavenumber, the "
+ "angular discretisation (phi in [0, 2*Pi] and theta in [0, Pi]) "
+ "of the far field sphere and the indices of the views containing the "
+ "complex-valued E and H fields. If `Normalize' is set, the far field "
+ "is normalized to 1. If `dB' is set, the far field is computed in dB. "
+ "If `NegativeTime' is set, E and H are assumed to have exp(-iwt) time "
+ "dependency; otherwise they are assume to have exp(+iwt) time "
+ "dependency.\n\n"
"Plugin(NearToFarField) creates one new view.";
}
@@ -46,11 +59,14 @@ StringXNumber *GMSH_NearToFarFieldPlugin::getOption(int iopt)
return &NearToFarFieldOptions_Number[iopt];
}
-double GMSH_NearToFarFieldPlugin::getFarField(std::vector<element*> allElems,
- std::vector<std::vector<double> > js,
- std::vector<std::vector<double> > ms,
- double k0, double r_far, double theta,
- double phi)
+// Compute field using e^{j\omega t} time dependency, following Jin in "Finite
+// Element Analysis of Antennas and Arrays", p. 176. This is not the usual `far
+// field', as it still contains the e^{ikr}/r factor.
+double GMSH_NearToFarFieldPlugin::getFarFieldJin(std::vector<element*> &allElems,
+ std::vector<std::vector<double> > &js,
+ std::vector<std::vector<double> > &ms,
+ double k0, double rFar, double theta,
+ double phi)
{
// theta in [0, pi] (elevation/polar angle)
// phi in [0, 2*pi] (azimuthal angle)
@@ -70,8 +86,8 @@ double GMSH_NearToFarFieldPlugin::getFarField(std::vector<element*> allElems,
N.resize(numSteps); Ns.resize(numSteps);
L.resize(numSteps); Ls.resize(numSteps);
for (int step = 0; step < numSteps; step++){
- N[step].resize(numComps); Ns[step].resize(numComps);
- L[step].resize(numComps); Ls[step].resize(numComps);
+ N[step].resize(numComps, 0.); Ns[step].resize(numComps);
+ L[step].resize(numComps, 0.); Ls[step].resize(numComps);
}
int i = 0 ;
@@ -123,9 +139,9 @@ double GMSH_NearToFarFieldPlugin::getFarField(std::vector<element*> allElems,
// E_r radial component is negligible in far field
double E_theta[2] ;
double E_phi[2] ;
- double k0_over_4pir = k0/(4*M_PI*r_far) ;
- double cos_k0r = cos(k0*r_far) ;
- double sin_k0r = sin(k0*r_far) ;
+ double k0_over_4pir = k0/(4*M_PI*rFar) ;
+ double cos_k0r = cos(k0*rFar) ;
+ double sin_k0r = sin(k0*rFar) ;
// Elevation component
E_theta[0] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * sin_k0r -
@@ -144,16 +160,77 @@ double GMSH_NearToFarFieldPlugin::getFarField(std::vector<element*> allElems,
return farF ;
}
+// Compute far field using e^{-i\omega t} time dependency, following Monk in
+// "Finite Element Methods for Maxwell's equations", p. 233
+double GMSH_NearToFarFieldPlugin::getFarFieldMonk(std::vector<element*> &allElems,
+ std::vector<std::vector<double> > &js,
+ std::vector<std::vector<double> > &ms,
+ double k0, double theta, double phi)
+{
+ double sTheta = sin(theta) ; double cTheta = cos(theta) ;
+ double sPhi = sin(phi) ; double cPhi = cos(phi) ;
+ double xHat[3] = {sTheta * cPhi, sTheta * sPhi, cTheta};
+ std::complex<double> I(0., 1.);
+ double Z0 = 120 * M_PI ; // free-space impedance
+
+ double integral_r[3] = {0., 0., 0.}, integral_i[3] = {0., 0., 0.};
+ int i = 0 ;
+ for(unsigned int ele = 0; ele < allElems.size(); ele++){
+ element* e = allElems[ele] ;
+ int numNodes = e->getNumNodes() ;
+ std::vector<double> integrand_r(numNodes * 3), integrand_i(numNodes * 3);
+ for(int nod = 0; nod < numNodes; nod++){
+ double y[3], xHat_dot_y;
+ e->getXYZ(nod, y[0], y[1], y[2]) ;
+ prosca(xHat, y, &xHat_dot_y) ;
+ double n_x_e_r[3] = {-ms[0][i], -ms[0][i + 1], -ms[0][i + 2]};
+ double n_x_e_i[3] = {-ms[1][i], -ms[1][i + 1], -ms[1][i + 2]};
+ double n_x_h_r[3] = {js[0][i], js[0][i + 1], js[0][i + 2]};
+ double n_x_h_i[3] = {js[1][i], js[1][i + 1], js[1][i + 2]};
+ double n_x_h_x_xHat_r[3], n_x_h_x_xHat_i[3];
+ prodve(n_x_h_r, xHat, n_x_h_x_xHat_r);
+ prodve(n_x_h_i, xHat, n_x_h_x_xHat_i);
+ for(int comp = 0; comp < 3; comp++){
+ std::complex<double> n_x_e(n_x_e_r[comp], n_x_e_i[comp]);
+ std::complex<double> n_x_h_x_xHat(n_x_h_x_xHat_r[comp], n_x_h_x_xHat_i[comp]);
+ // Warning: Z0 == 1 in Monk
+ std::complex<double> integrand = (n_x_e + Z0 * n_x_h_x_xHat) *
+ (cos(-k0 * xHat_dot_y) + I * sin(-k0 * xHat_dot_y));
+ integrand_r[3 * nod + comp] = integrand.real();
+ integrand_i[3 * nod + comp] = integrand.imag();
+ }
+ i += 3;
+ }
+ for(int comp = 0; comp < 3; comp++){
+ integral_r[comp] += e->integrate(&integrand_r[comp], 3);
+ integral_i[comp] += e->integrate(&integrand_i[comp], 3);
+ }
+ }
+
+ double xHat_x_integral_r[3], xHat_x_integral_i[3];
+ prodve(xHat, integral_r, xHat_x_integral_r);
+ prodve(xHat, integral_i, xHat_x_integral_i);
+ std::complex<double> coef = I * k0 / 4. / M_PI;
+ std::complex<double> einf[3] = {coef * (xHat_x_integral_r[0] + I * xHat_x_integral_i[0]),
+ coef * (xHat_x_integral_r[1] + I * xHat_x_integral_i[1]),
+ coef * (xHat_x_integral_r[2] + I * xHat_x_integral_i[2])};
+ return (norm(einf[0]) + norm(einf[1]) + norm(einf[2]));
+}
+
PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
{
double _k0 = (double)NearToFarFieldOptions_Number[0].def;
- double _r_far = (double)NearToFarFieldOptions_Number[1].def;
- int _NbPhi = (int)NearToFarFieldOptions_Number[2].def;
- int _NbThe = (int)NearToFarFieldOptions_Number[3].def;
- int _eView = (int)NearToFarFieldOptions_Number[4].def;
- int _hView = (int)NearToFarFieldOptions_Number[5].def;
- bool _normalize = (bool)NearToFarFieldOptions_Number[6].def;
- bool _dB = (bool)NearToFarFieldOptions_Number[7].def;
+ double _phiStart = (double)NearToFarFieldOptions_Number[1].def;
+ double _phiEnd = (double)NearToFarFieldOptions_Number[2].def;
+ int _NbPhi = (int)NearToFarFieldOptions_Number[3].def;
+ double _thetaStart = (double)NearToFarFieldOptions_Number[4].def;
+ double _thetaEnd = (double)NearToFarFieldOptions_Number[5].def;
+ int _NbThe = (int)NearToFarFieldOptions_Number[6].def;
+ int _eView = (int)NearToFarFieldOptions_Number[7].def;
+ int _hView = (int)NearToFarFieldOptions_Number[8].def;
+ bool _normalize = (bool)NearToFarFieldOptions_Number[9].def;
+ bool _dB = (bool)NearToFarFieldOptions_Number[10].def;
+ int _negativeTime = (int)NearToFarFieldOptions_Number[11].def;
PView *ve = getView(_eView, v);
if(!ve){
@@ -171,34 +248,33 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
if(eData->getNumEntities() != hData->getNumEntities() ||
eData->getNumElements() != hData->getNumElements() ||
eData->getNumTimeSteps()!= hData->getNumTimeSteps()){
- Msg::Error("Incompatible views for e-field and h-field");
+ Msg::Error("Incompatible views for E field and H field");
return v;
}
if(eData->getNumTimeSteps() != 2 || hData->getNumTimeSteps() != 2){
- Msg::Error("Invalid number of steps for EView or HView, fields must be complex");
+ Msg::Error("Invalid number of steps for E or H fields (must be complex)");
return v;
}
// center of the Far Field sphere
- double x0 = eData->getBoundingBox().center().x();
- double y0 = eData->getBoundingBox().center().y();
- double z0 = eData->getBoundingBox().center().z();
+ SBoundingBox3d bbox = eData->getBoundingBox();
+ double x0 = bbox.center().x();
+ double y0 = bbox.center().y();
+ double z0 = bbox.center().z();
+ double lc = norm(SVector3(bbox.max(), bbox.min()));
+
if(x0 != hData->getBoundingBox().center().x() ||
y0 != hData->getBoundingBox().center().y() ||
z0 != hData->getBoundingBox().center().z()){
- Msg::Error("EView %i and HView %i must be given on the same grid", _eView, _hView);
+ Msg::Error("E and H fields must be given on the same grid");
return v;
}
// compute surface currents on all input elements
std::vector<element*> allElems ;
- std::vector<std::vector<double> > js ;
- std::vector<std::vector<double> > ms ;
-
- int numSteps = eData->getNumTimeSteps() ;
- js.resize(numSteps);
- ms.resize(numSteps);
+ std::vector<std::vector<double> > js(2) ;
+ std::vector<std::vector<double> > ms(2) ;
for(int ent = 0; ent < eData->getNumEntities(0); ent++){
for(int ele = 0; ele < eData->getNumElements(0, ent); ele++){
@@ -222,7 +298,7 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
else
normal2points(x[0], y[0], z[0], x[1], y[1], z[1], n);
- for(int step = 0; step < numSteps; step++){
+ for(int step = 0; step < 2; step++){
for(int nod = 0; nod < numNodes; nod++){
double h[3], e[3];
for(int comp = 0; comp < numComp; comp++){
@@ -249,50 +325,49 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
PView *vf = new PView();
PViewDataList *dataFar = getDataList(vf);
- double phi, dPhi = 2*M_PI/_NbPhi ;
- double theta, dTheta = M_PI/_NbThe ;
- double ffmax = 0.0 ;
-
- std::vector<std::vector<double> > allPhi ;
- std::vector<std::vector<double> > allThe ;
- std::vector<std::vector<double> > farF ;
-
- allPhi.resize(_NbPhi+1);
- allThe.resize(_NbPhi+1);
- farF.resize(_NbPhi+1);
- for (int i = 0; i<=_NbPhi; i++){
- allPhi[i].resize(_NbThe+1);
- allThe[i].resize(_NbThe+1);
- farF[i].resize(_NbThe+1);
+ std::vector<std::vector<double> > allPhi(_NbPhi + 1);
+ std::vector<std::vector<double> > allThe(_NbPhi + 1);
+ std::vector<std::vector<double> > farF(_NbPhi + 1);
+ for (int i = 0; i <= _NbPhi; i++){
+ allPhi[i].resize(_NbThe + 1);
+ allThe[i].resize(_NbThe + 1);
+ farF[i].resize(_NbThe + 1);
}
+ double dPhi = (_phiEnd - _phiStart) / _NbPhi ;
+ double dTheta = (_thetaEnd - _thetaStart) / _NbThe ;
+ double ffmax = 0.0 ;
+ Msg::ResetProgressMeter();
for (int i = 0; i <= _NbPhi; i++){
- phi = i*dPhi ;
+ double phi = _phiStart + i * dPhi ;
for (int j = 0; j <= _NbThe; j++){
- theta = j * dTheta ;
+ double theta = _thetaStart + j * dTheta ;
allPhi[i][j] = phi ;
allThe[i][j] = theta ;
- farF[i][j] = getFarField(allElems, js, ms, _k0, _r_far, theta, phi) ;
+ if(_negativeTime)
+ farF[i][j] = getFarFieldMonk(allElems, js, ms, _k0, theta, phi);
+ else
+ farF[i][j] = getFarFieldJin(allElems, js, ms, _k0, 10 * lc, theta, phi);
ffmax = (ffmax < farF[i][j]) ? farF[i][j] : ffmax ;
}
+ Msg::ProgressMeter(i, _NbPhi, "Computing far field");
}
if(_normalize){
for (int i = 0; i <= _NbPhi; i++)
for (int j = 0; j <= _NbThe; j++)
- if(ffmax!=0.0)
+ if(ffmax != 0.0)
farF[i][j] /= ffmax ;
else
Msg::Warning("Far field pattern not normalized, max value = %g", ffmax);
}
+ for(unsigned int i = 0; i < allElems.size(); i++)
+ delete allElems[i];
+
// construct sphere for visualization, centered at center of bb and with
// radius relative to the bb size
- double r_bb[3] = {eData->getBoundingBox().max().x()-eData->getBoundingBox().min().x(),
- eData->getBoundingBox().max().y()-eData->getBoundingBox().min().y(),
- eData->getBoundingBox().max().z()-eData->getBoundingBox().min().z()};
- double r_sph = norm3(r_bb) ;
- r_sph = (r_sph) ? r_sph/2 : 1./2. ; // radious of sphere for visu
+ double r_sph = lc ? lc / 2. : 1; // radius of sphere for visu
for (int i = 0; i < _NbPhi; i++){
for (int j = 0; j < _NbThe; j++){
diff --git a/Plugin/NearToFarField.h b/Plugin/NearToFarField.h
index 618982623d..309608d0ad 100644
--- a/Plugin/NearToFarField.h
+++ b/Plugin/NearToFarField.h
@@ -28,13 +28,20 @@ class GMSH_NearToFarFieldPlugin : public GMSH_PostPlugin
return "Compute Far Field pattern from Near Field on a surface";
}
std::string getHelp() const;
+ virtual std::string getAuthor() const { return "R. Sabariego, C. Geuzaine"; }
int getNbOptions() const;
StringXNumber* getOption(int iopt);
PView *execute(PView *);
- static double getFarField(std::vector<element*> allElems,
- std::vector<std::vector<double> >js, std::vector<std::vector<double> >ms,
- double k0, double r_far, double theta, double phi);
- };
+ double getFarFieldJin(std::vector<element*> &allElems,
+ std::vector<std::vector<double> > &js,
+ std::vector<std::vector<double> > &ms,
+ double k0, double r_far, double theta, double phi);
+
+ double getFarFieldMonk(std::vector<element*> &allElems,
+ std::vector<std::vector<double> > &js,
+ std::vector<std::vector<double> > &ms,
+ double k0, double theta, double phi);
+};
#endif
diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h
index 8e056f2d20..1ebb8d0837 100644
--- a/Post/shapeFunctions.h
+++ b/Post/shapeFunctions.h
@@ -395,17 +395,22 @@ public:
default: start = end = 0; break;
}
}
- inline int getNumGaussPoints(){ return 3; }
+ inline int getNumGaussPoints(){ return /* 3 */ 1; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u3[3] = {0.16666666666666,0.66666666666666,0.16666666666666};
- double v3[3] = {0.16666666666666,0.16666666666666,0.66666666666666};
- 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];
+ /*
+ 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];
+ */
+ if(num < 0 || num > 0) return;
+ u = v = 0.333333333333333; w = 0.;
+ weight = 0.5;
}
void getShapeFunction(int num, double u, double v, double w, double &s)
{
@@ -491,9 +496,9 @@ public:
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u4[4] = {0.577350269189,-0.577350269189,0.577350269189,-0.577350269189};
- double v4[4] = {0.577350269189,0.577350269189,-0.577350269189,-0.577350269189};
- double p4[4] = {1.,1.,1.,1.};
+ 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];
@@ -573,10 +578,10 @@ public:
inline int getNumGaussPoints(){ return 4; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u4[4] = {0.138196601125,0.138196601125,0.138196601125,0.585410196625};
- double v4[4] = {0.138196601125,0.138196601125,0.585410196625,0.138196601125};
- double w4[4] = {0.138196601125,0.585410196625,0.138196601125,0.138196601125};
- double p4[4] = {0.0416666666667,0.0416666666667,0.0416666666667,0.0416666666667};
+ 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];
@@ -670,14 +675,14 @@ public:
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u6[6] = { 0.40824826, 0.40824826, -0.40824826,
- -0.40824826, -0.81649658, 0.81649658};
- double v6[6] = { 0.70710678, -0.70710678, 0.70710678,
- -0.70710678, 0., 0.};
- double w6[6] = {-0.57735027, -0.57735027, 0.57735027,
- 0.57735027, -0.57735027, 0.57735027};
- double p6[6] = { 1.3333333333, 1.3333333333, 1.3333333333,
- 1.3333333333, 1.3333333333, 1.3333333333};
+ 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];
@@ -772,14 +777,14 @@ public:
inline int getNumGaussPoints(){ return 6; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u6[6] = {0.166666666666666, 0.333333333333333, 0.166666666666666,
- 0.166666666666666, 0.333333333333333, 0.166666666666666};
- double v6[6] = {0.166666666666666, 0.166666666666666, 0.333333333333333,
- 0.166666666666666, 0.166666666666666, 0.333333333333333};
- double w6[6] = {-0.577350269189, -0.577350269189, -0.577350269189,
- 0.577350269189, 0.577350269189, 0.577350269189};
- double p6[6] = {0.166666666666666,0.166666666666666,0.166666666666666,
- 0.166666666666666,0.166666666666666,0.166666666666666,};
+ 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];
@@ -864,22 +869,22 @@ public:
inline int getNumGaussPoints(){ return 8; }
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
{
- double u8[8] = {0.2631840555694285,-0.2631840555694285,
- 0.2631840555694285,-0.2631840555694285,
- 0.5066163033492386,-0.5066163033492386,
- 0.5066163033492386,-0.5066163033492386};
- double v8[8] = {0.2631840555694285,0.2631840555694285,
- -0.2631840555694285,-0.2631840555694285,
- 0.5066163033492386,0.5066163033492386,
- -0.5066163033492386,-0.5066163033492386};
- double w8[8] = {0.544151844011225,0.544151844011225,
- 0.544151844011225,0.544151844011225,
- 0.122514822655441,0.122514822655441,
- 0.122514822655441,0.122514822655441};
- double p8[8] = {0.100785882079825,0.100785882079825,
- 0.100785882079825,0.100785882079825,
- 0.232547451253508,0.232547451253508,
- 0.232547451253508,0.232547451253508};
+ static double u8[8] = {0.2631840555694285,-0.2631840555694285,
+ 0.2631840555694285,-0.2631840555694285,
+ 0.5066163033492386,-0.5066163033492386,
+ 0.5066163033492386,-0.5066163033492386};
+ static double v8[8] = {0.2631840555694285,0.2631840555694285,
+ -0.2631840555694285,-0.2631840555694285,
+ 0.5066163033492386,0.5066163033492386,
+ -0.5066163033492386,-0.5066163033492386};
+ 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.100785882079825,0.100785882079825,
+ 0.100785882079825,0.100785882079825,
+ 0.232547451253508,0.232547451253508,
+ 0.232547451253508,0.232547451253508};
if(num < 0 || num > 7) return;
u = u8[num];
v = v8[num];
--
GitLab