Commit 44f4bc8c by Christophe Geuzaine

cleanup

parent 3f1a8068
......@@ -484,9 +484,9 @@ bool OCCFace::containsParam(const SPoint2 &pt)
SPoint2 gp2 = stl_vertices[stl_triangles[i + 1]];
SPoint2 gp3 = stl_vertices[stl_triangles[i + 2]];
double s1 = robustPredicates::orient2d (gp1,gp2,mine);
double s2 = robustPredicates::orient2d (gp2,gp3,mine);
double s3 = robustPredicates::orient2d (gp3,gp1,mine);
double s1 = robustPredicates::orient2d(gp1, gp2, mine);
double s2 = robustPredicates::orient2d(gp2, gp3, mine);
double s3 = robustPredicates::orient2d(gp3, gp1, mine);
/*
fprintf(F,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
gp1.x(),gp1.y(),0.0,
......
......@@ -12,20 +12,18 @@
class bezierBasis;
class GradientBasis {
public:
public:
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
fullMatrix<double> gradShapeIdealMatX, gradShapeIdealMatY, gradShapeIdealMatZ;
private:
private:
const FuncSpaceData _data;
public:
public:
GradientBasis(FuncSpaceData);
int getNumSamplingPoints() const {return gradShapeMatX.size1();}
int getNumMapNodes() const {return gradShapeMatX.size2();}
const bezierBasis* getBezier() const;
void getGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
......@@ -40,13 +38,17 @@ public:
fullMatrix<double> &dxyzdXYZ) const;
void mapFromIdealElement(fullMatrix<double> &dxyzdX,
fullMatrix<double> &dxyzdY,
fullMatrix<double> &dxyzdZ) const {
GradientBasis::mapFromIdealElement(_data.elementType(), dxyzdX, dxyzdY, dxyzdZ);
fullMatrix<double> &dxyzdZ) const
{
GradientBasis::mapFromIdealElement(_data.elementType(),
dxyzdX, dxyzdY, dxyzdZ);
}
void mapFromIdealElement(fullVector<double> &dxyzdX,
fullVector<double> &dxyzdY,
fullVector<double> &dxyzdZ) const {
GradientBasis::mapFromIdealElement(_data.elementType(), dxyzdX, dxyzdY, dxyzdZ);
fullVector<double> &dxyzdZ) const
{
GradientBasis::mapFromIdealElement(_data.elementType(),
dxyzdX, dxyzdY, dxyzdZ);
}
static void mapFromIdealElement(int type,
fullMatrix<double> &gSMatX,
......@@ -57,27 +59,25 @@ public:
fullVector<double> &gSVecY,
fullVector<double> &gSVecZ);
static void mapFromIdealElement(int type, double jac[3][3]);
void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
};
class JacobianBasis {
private:
private:
const GradientBasis *_gradBasis;
const FuncSpaceData _data;
const int _dim;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBaryX, primGradShapeBaryY, primGradShapeBaryZ;
fullVector<double> primIdealGradShapeBaryX, primIdealGradShapeBaryY,
primIdealGradShapeBaryZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
fullVector<double> primIdealGradShapeBaryX, primIdealGradShapeBaryY;
fullVector<double> primIdealGradShapeBaryZ;
// Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
fullMatrix<double> matrixPrimJac2Jac;
int numJacNodes, numPrimJacNodes;
int numMapNodes, numPrimMapNodes;
int numJacNodesFast;
public:
public:
JacobianBasis(FuncSpaceData);
// Get methods
......@@ -97,7 +97,8 @@ public:
double getPrimJac3D(const fullMatrix<double> &nodesXYZ, bool ideal=false) const;
inline void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const {
fullMatrix<double> &JDJ) const
{
getSignedJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
......@@ -112,7 +113,8 @@ public:
}
inline void getSignedIdealJacAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const {
fullMatrix<double> &JDJ) const
{
getSignedJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
_gradBasis->gradShapeIdealMatY,
_gradBasis->gradShapeIdealMatZ,
......@@ -124,7 +126,8 @@ public:
fullMatrix<double> &gradLambdaJ) const;
inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
......@@ -134,7 +137,8 @@ public:
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
fullMatrix<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
......@@ -142,7 +146,8 @@ public:
}
inline void getSignedIdealJacobian(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
_gradBasis->gradShapeIdealMatY,
_gradBasis->gradShapeIdealMatZ,
......@@ -152,14 +157,16 @@ public:
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
fullMatrix<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
_gradBasis->gradShapeIdealMatY,
_gradBasis->gradShapeIdealMatZ,
nodesX, nodesY, nodesZ, true, false, jacobian, normals);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian) const {
fullVector<double> &jacobian) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
......@@ -168,7 +175,8 @@ public:
inline void getScaledJacobian(const fullMatrix<double> &nodesX,
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
fullMatrix<double> &jacobian) const {
fullMatrix<double> &jacobian) const
{
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
......@@ -176,14 +184,16 @@ public:
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
gradShapeMatYFast, gradShapeMatZFast,
nodesXYZ, false, false, jacobian, normals);
}
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian,
const fullMatrix<double> *normals = NULL) const {
const fullMatrix<double> *normals = NULL) const
{
getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
gradShapeMatYFast, gradShapeMatZFast,
nodesXYZ, false, true, jacobian, normals);
......@@ -191,7 +201,9 @@ public:
void lag2Bez(const fullVector<double> &lag, fullVector<double> &bez) const;
void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
inline void primJac2Jac(const fullVector<double> &primJac, fullVector<double> &jac) const {
inline void primJac2Jac(const fullVector<double> &primJac,
fullVector<double> &jac) const
{
matrixPrimJac2Jac.mult(primJac,jac);
}
......@@ -199,12 +211,10 @@ public:
void interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
fullMatrix<double> &result, bool areBezier = false) const;
static int jacobianOrder(int tag);
static int jacobianOrder(int parentType, int order);
static FuncSpaceData jacobianMatrixSpace(int type, int order);
private :
void getJacobianGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
......@@ -224,7 +234,6 @@ public:
bool idealNorm, bool scaling,
fullMatrix<double> &jacobian,
const fullMatrix<double> *normals) const;
void getSignedJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
......
......@@ -86,27 +86,28 @@ class bezierBasis {
class bezierBasisRaiser {
// Let f, g, h be three function whose Bezier coefficients are given.
// This class allows to compute the Bezier coefficients of f*g and f*g*h.
private :
private :
class _Data {
friend class bezierBasisRaiser;
private:
private:
int i, j, k;
double val;
public:
public:
_Data(double vv, int ii, int jj = -1, int kk = -1) :
i(ii), j(jj), k(kk), val(vv) {}
};
std::vector<std::vector<_Data> > _raiser2, _raiser3;
const bezierBasis *_bfs;
public:
bezierBasisRaiser(const bezierBasis *bezier) : _bfs(bezier) {
public:
bezierBasisRaiser(const bezierBasis *bezier) : _bfs(bezier)
{
_fillRaiserData();
};
const bezierBasis* getRaisedBezierBasis(int multipliedBy2or3);
// const bezierBasis* getRaisedBezierBasis(int raised) const;
// const bezierBasis* getRaisedBezierBasis(int raised) const;
void computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffB,
......@@ -122,13 +123,9 @@ public:
const fullMatrix<double> &coeffB,
const fullMatrix<double> &coeffC,
fullMatrix<double> &coeffCubic);
private:
private:
void _fillRaiserData();
void _fillRaiserDataPyr();
};
#endif
......@@ -18,12 +18,13 @@ class nodalBasis {
nodalBasis() {}
nodalBasis(int tag);
virtual ~nodalBasis() {}
virtual int getNumShapeFunctions() const = 0;
void getReferenceNodes(fullMatrix<double> &nodes) const {
void getReferenceNodes(fullMatrix<double> &nodes) const
{
nodes = points;
}
const fullMatrix<double>& getReferenceNodes() const {
const fullMatrix<double>& getReferenceNodes() const
{
return points;
}
void getReferenceNodesForBezier(fullMatrix<double> &nodes) const;
......@@ -33,17 +34,20 @@ class nodalBasis {
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const = 0;
virtual void df(double u, double v, double w, double grads[][3]) const = 0;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const = 0;
virtual void ddf(double u, double v, double w, double grads[][3][3]) const {Msg::Fatal("Not implemented");}
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {Msg::Fatal("Not implemented");}
virtual void ddf(double u, double v, double w, double grads[][3][3]) const {
Msg::Error("ddf not implemented for this basis");
}
virtual void dddf(double u, double v, double w, double grads[][3][3][3]) const {
Msg::Error("dddf not implemented for this basis");
}
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
// nodes of the element that rotates the element so that the
// considered face becomes the first one in the right orientation;
// For element, like prisms that have different kind of faces,
// fullCLosure[i] rotates the element so that the considered face
// becomes the closureRef[i]-th face (the first tringle or the first
// quad face)
// closures is the list of the nodes of each face, in the order of the
// polynomialBasis of the face; fullClosures is mapping of the nodes of the
// element that rotates the element so that the considered face becomes the
// first one in the right orientation; For element, like prisms that have
// different kind of faces, fullCLosure[i] rotates the element so that the
// considered face becomes the closureRef[i]-th face (the first triangle or
// the first quad face)
class closure : public std::vector<int> {
public:
int type;
......@@ -52,11 +56,17 @@ class nodalBasis {
clCont closures, fullClosures;
std::vector<int> closureRef;
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
// for a given face/edge, with both a sign and a rotation, give an ordered
// list of nodes on this face/edge
virtual int getClosureType(int id) const { return closures[id].type; }
virtual const std::vector<int> &getClosure(int id) const { return closures[id]; }
virtual const std::vector<int> &getFullClosure(int id) const { return fullClosures[id]; }
virtual const std::vector<int> &getClosure(int id) const
{
return closures[id];
}
virtual const std::vector<int> &getFullClosure(int id) const
{
return fullClosures[id];
}
inline int getClosureId(int iFace, int iSign = 1, int iRot = 0) const;
inline void breakClosureId(int i, int &iFace, int &iSign, int &iRot) const;
};
......
......@@ -8,13 +8,10 @@
#include "fullMatrix.h"
#include "FuncSpaceData.h"
/*
* Functions to generate point distributions on
* the references elements, for all orders.
* &
* Functions generating exponents of Pascal monomials
* in the same order than Gmsh Points.
*/
// Functions to generate point distributions on the references elements, for all
// orders and functions generating exponents of Pascal monomials in the same
// order than Gmsh Points.
// Points
......@@ -30,8 +27,8 @@ fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip = false
fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip = false);
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip = false);
fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, bool serendip = false);
fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk,
bool serendip = false);
// Monomial exponents
......@@ -40,23 +37,29 @@ void gmshGenerateMonomials(FuncSpaceData, fullMatrix<double> &);
fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip = false);
fullMatrix<double> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
fullMatrix<double> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsQuadrangle(int order, bool
forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsQuadSerendipity(int order);
fullMatrix<double> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false);
fullMatrix<double> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsHexahedron(int order,
bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order);
fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order);
// Generate monomials of pyramidal nodal space {X^i Y^j Z^k | i,j <= k, k <= 'order'},
fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
// Generate monomials of pyramidal nodal space {X^i Y^j Z^k | i,j <= k, k <=
// 'order'},
fullMatrix<double> gmshGenerateMonomialsPyramid(int order,
bool forSerendipPoints = false);
fullMatrix<double> gmshGenerateMonomialsPyramidSerendipity(int order);
// If 'serendip' == true, generate monomials of serendipity pyramid at order 'nk',
// else if 'pyr' == true, generate monomials of space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
// else if 'pyr' == false, generate monomials of space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij, int nk, bool forSerendipPoints = false);
// If 'serendip' == true,
// generate monomials of serendipity pyramid at order 'nk',
// else if 'pyr' == true,
// generate monomials of space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
// else if 'pyr' == false,
// generate monomials of space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij, int nk,
bool forSerendipPoints = false);
#endif
......@@ -13,7 +13,6 @@
#include "nodalBasis.h"
#include <iostream>
inline double pow_int(const double &a, const int &n)
{
switch (n) {
......@@ -73,8 +72,8 @@ inline double pow_int(const double &a, const double &d)
class polynomialBasis : public nodalBasis
{
public:
// for now the only implemented polynomial basis are nodal poly
// basis, we use the type of the corresponding gmsh element as type
// for now the only implemented polynomial basis are nodal poly basis, we use
// the type of the corresponding gmsh element as type
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
......@@ -82,7 +81,7 @@ class polynomialBasis : public nodalBasis
polynomialBasis(int tag);
~polynomialBasis();
virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
virtual inline int getNumShapeFunctions() const { return coefficients.size1(); }
virtual void f(double u, double v, double w, double *sf) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
......@@ -94,5 +93,4 @@ class polynomialBasis : public nodalBasis
void evaluateMonomials(double u, double v, double w, double p[]) const;
};
#endif
......@@ -11,10 +11,11 @@
pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag), bergot(0)
{
if (serendip && order > 2) {
Msg::Warning("Serendipity pyramid for order %i not yet implemented",order);
Msg::Warning("Serendipity pyramid for order %i not yet implemented",
order);
return;
}
bergot = new BergotBasis(order,serendip);
int num_points = points.size1();
......@@ -41,7 +42,7 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag), bergot(0)
monomials(idx,0) = i;
monomials(idx,1) = j;
monomials(idx,2) = k;
for (int l=0;l<num_points;l++) {
double oneMinW = std::max(1e-14,1-points(l,2));
VDM(idx,l) = std::pow(points(l,0),i);
......
......@@ -18,19 +18,16 @@ class pyramidalBasis: public nodalBasis
fullMatrix<double> bergotCoefficients;
public:
fullMatrix<double> coefficients;
fullMatrix<double> monomials;
public:
pyramidalBasis(int tag);
~pyramidalBasis();
virtual void f(double u, double v, double w, double *val) const;
virtual void f(const fullMatrix<double> &coord, fullMatrix<double> &sf) const;
virtual void df(double u, double v, double w, double grads[][3]) const;
virtual void df(const fullMatrix<double> &coord, fullMatrix<double> &dfm) const;
virtual int getNumShapeFunctions() const;
};
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment