From 01fef23fcc80359a95f598abd32ee1d5e35ddba2 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 29 Aug 2010 17:48:55 +0000
Subject: [PATCH] FIXME: please remove dependencies to Geo/ code in Numeric/ !
---
Geo/partitionEdge.h | 5 ++---
Geo/partitionFace.h | 5 ++---
Geo/partitionVertex.h | 4 ++--
Numeric/DivideAndConquer.cpp | 2 ++
Numeric/DivideAndConquer.h | 2 ++
Numeric/Numeric.cpp | 14 +++-----------
Numeric/Numeric.h | 26 ++++++++++++++------------
Numeric/fullMatrix.cpp | 3 ++-
Numeric/fullMatrix.h | 12 ++++++++----
Numeric/polynomialBasis.cpp | 34 +++++++++++++++++++++-------------
Numeric/polynomialBasis.h | 14 +++++++++-----
Numeric/simpleFunction.h | 18 ++++++++++--------
12 files changed, 77 insertions(+), 62 deletions(-)
diff --git a/Geo/partitionEdge.h b/Geo/partitionEdge.h
index e9324ac82b..8888e6aa0c 100644
--- a/Geo/partitionEdge.h
+++ b/Geo/partitionEdge.h
@@ -16,15 +16,14 @@ class partitionEdge : public discreteEdge {
public:
partitionEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1,
std::vector<int> &partitions)
- : discreteEdge(model,num,_v0,_v1),_partitions(partitions)
+ : discreteEdge(model, num, _v0, _v1), _partitions(partitions)
{
- std::sort(_partitions.begin(),_partitions.end());
+ std::sort(_partitions.begin(), _partitions.end());
}
virtual ~partitionEdge() {}
virtual GeomType geomType() const { return PartitionCurve; }
};
-
struct Less_partitionEdge :
public std::binary_function<partitionEdge*, partitionEdge*, bool> {
bool operator()(const partitionEdge* e1, const partitionEdge* e2) const
diff --git a/Geo/partitionFace.h b/Geo/partitionFace.h
index 85e89049ff..d2971b0a2c 100644
--- a/Geo/partitionFace.h
+++ b/Geo/partitionFace.h
@@ -14,15 +14,14 @@ class partitionFace : public discreteFace {
std::vector<int> _partitions;
public:
partitionFace(GModel *model, int num, std::vector<int> &partitions)
- : discreteFace(model,num),_partitions(partitions)
+ : discreteFace(model, num), _partitions(partitions)
{
- std::sort(_partitions.begin(),_partitions.end());
+ std::sort(_partitions.begin(), _partitions.end());
}
virtual ~partitionFace() {}
virtual GeomType geomType() const { return PartitionSurface; }
};
-
struct Less_partitionFace :
public std::binary_function<partitionFace*, partitionFace*, bool> {
bool operator()(const partitionFace* e1, const partitionFace* e2) const
diff --git a/Geo/partitionVertex.h b/Geo/partitionVertex.h
index 3f25db98ec..7336affb25 100644
--- a/Geo/partitionVertex.h
+++ b/Geo/partitionVertex.h
@@ -14,9 +14,9 @@ class partitionVertex : public discreteVertex {
std::vector<int> _partitions;
public:
partitionVertex(GModel *model, int num, std::vector<int> &partitions)
- : discreteVertex(model,num),_partitions(partitions)
+ : discreteVertex(model, num), _partitions(partitions)
{
- std::sort(_partitions.begin(),_partitions.end());
+ std::sort(_partitions.begin(), _partitions.end());
}
virtual ~partitionVertex() {}
virtual GeomType geomType() const { return PartitionVertex; }
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 76dc26fd55..38c4607cdc 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -19,6 +19,8 @@
#include "DivideAndConquer.h"
#include "Numeric.h"
#include "robustPredicates.h"
+
+// FIXME: should not introduce dependencies to Geo/ code in Numeric
#include "GPoint.h"
#include "GFace.h"
#include "MLine.h"
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 223426b60d..d3008f45ba 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -9,6 +9,8 @@
#include <vector>
#include <algorithm>
#include "fullMatrix.h"
+
+// FIXME: should not introduce dependencies to Geo/ code in Numeric
#include "SPoint2.h"
#include "simpleFunction.h"
#include "Octree.h"
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 7ff6329166..14fd1cc780 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -16,7 +16,6 @@ double myatan2(double a, double b)
return atan2(a, b);
}
-
double myasin(double a)
{
if(a <= -1.)
@@ -80,7 +79,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
if(norm == 0.0 || fabs(det) / norm < 1.e-12) {
if(norm)
Msg::Debug("Assuming 2x2 matrix is singular (det/norm == %.16g)",
- fabs(det) / norm);
+ fabs(det) / norm);
res[0] = res[1] = 0.0;
return 0;
}
@@ -162,7 +161,7 @@ int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
if(norm == 0.0 || fabs(*det) / norm < 1.e-12) {
if(norm)
Msg::Debug("Assuming 3x3 matrix is singular (det/norm == %.16g)",
- fabs(*det) / norm);
+ fabs(*det) / norm);
res[0] = res[1] = res[2] = 0.0;
return 0;
}
@@ -290,12 +289,6 @@ double triangle_area2d(double p0[2], double p1[2], double p2[2])
return 0.5 * sqrt(c*c);
}
-double triangle_polar_inertia(double p0[2], double p1[2], double p2[2])
-{
- throw;
-}
-
-
void circumCenterXY(double *p1, double *p2, double *p3, double *res)
{
double d, a1, a2, a3;
@@ -309,7 +302,7 @@ void circumCenterXY(double *p1, double *p2, double *p3, double *res)
d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
if(d == 0.0) {
- // Msg::Warning("Colinear points in circum circle computation");
+ // Msg::Warning("Colinear points in circum circle computation");
res[0] = res[1] = -99999.;
return ;
}
@@ -321,7 +314,6 @@ void circumCenterXY(double *p1, double *p2, double *p3, double *res)
res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
}
-
void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv)
{
double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 8748b0da5e..33643dd618 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -63,7 +63,6 @@ double angle_02pi(double A3);
double angle_plan(double V[3], double P1[3], double P2[3], double n[3]);
double triangle_area(double p0[3], double p1[3], double p2[3]);
double triangle_area2d(double p0[2], double p1[2], double p2[2]);
-double triangle_polar_inertia(double p0[2], double p1[2], double p2[2]);
void circumCenterXY(double *p1, double *p2, double *p3, double *res);
void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv=0);
// operate a transformation on the 4 points of a Quad in 3D, to have an equivalent Quad in 2D
@@ -90,21 +89,24 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
const SPoint3 &p1,
const SPoint3 &p2,
const SPoint3 &p3);
-void signedDistancePointLine(const SPoint3 &p1,const SPoint3 &p2,const SPoint3 &p, double &distance, SPoint3 &closePt);
+void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p,
+ double &distance, SPoint3 &closePt);
void signedDistancesPointsLine (std::vector<double>&distances,
std::vector<SPoint3>&closePts,
const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2);
+ const SPoint3 &p1, const SPoint3 &p2);
-void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &closePt,const SPoint3 &p1,const SPoint3 &p2,double* xp,double* yp,double* otherp,double* x,double* y,double* other);
-int computeDistanceRatio(const double &y, const double &yp,const double &x,const double &xp, double *distance, const double &r1, const double &r2);
+void changeReferential(const int direction, const SPoint3 &p, const SPoint3 &closePt,
+ const SPoint3 &p1, const SPoint3 &p2, double *xp, double*yp,
+ double *otherp, double *x, double *y, double *other);
+int computeDistanceRatio(const double &y, const double &yp, const double &x,
+ const double &xp, double *distance, const double &r1,
+ const double &r2);
void signedDistancesPointsEllipseLine (std::vector<double>&distances,
- std::vector<double>&distancesE,
- std::vector<int>&isInYarn,
- std::vector<SPoint3>&closePts,
- const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2);
+ std::vector<double>&distancesE,
+ std::vector<int>&isInYarn,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1, const SPoint3 &p2);
#endif
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index bfc871fa1e..ea9d2a6661 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -329,7 +329,8 @@ void fullMatrix<double>::registerBindings(binding *b)
cm->setDescription("Sets the (i,j) entry of the matrix to v");
cm = cb->addMethod("resize", &fullMatrix<double>::resize);
cm->setArgNames("nRows","nColumns","reset",NULL);
- cm->setDescription("Change the size of the fullMatrix (and re-alloc if needed), values are set to zero if reset is true");
+ cm->setDescription("Change the size of the fullMatrix (and re-alloc if needed), "
+ "values are set to zero if reset is true");
cm = cb->addMethod("gemm", &fullMatrix<double>::gemm);
cm->setArgNames("A","B","alpha","beta",NULL);
cm->setDescription("this = beta*this + alpha * (A.B)");
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 250b189f92..fecedd9216 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -131,14 +131,16 @@ class fullMatrix
inline scalar get(int r, int c) const {
#ifdef _DEBUG
if (r >= _r || r < 0 || c >= _c || c < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)", r, c, _r, _c);
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ r, c, _r, _c);
#endif
return (*this)(r, c);
}
inline void set(int r, int c, scalar v){
#ifdef _DEBUG
if (r >= _r || r < 0 || c >= _c || c < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)", r, c, _r, _c);
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ r, c, _r, _c);
#endif
(*this)(r, c) = v;
}
@@ -260,7 +262,8 @@ class fullMatrix
{
#ifdef _DEBUG
if (i >= _r || i < 0 || j >= _c || j < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)", i, j, _r, _c);
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ i, j, _r, _c);
#endif
return _data[i + _r * j];
}
@@ -268,7 +271,8 @@ class fullMatrix
{
#ifdef _DEBUG
if (i >= _r || i < 0 || j >= _c || j < 0)
- Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)", i, j, _r, _c);
+ Msg::Fatal("invalid index to access fullMatrix : %i %i (size = %i %i)",
+ i, j, _r, _c);
#endif
return _data[i + _r * j];
}
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 6fc0e59235..5b91118476 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -93,6 +93,7 @@ static fullMatrix<double> generatePascalQuad(int order)
}
return monomials;
}
+
/*
00 10 20 30 40 站ｯ
01 11 21 31 41 站ｯ
@@ -795,7 +796,8 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
fi = fi + 3; ti = ti + 3;
for (int l = 0; l < orderint - 1; l++){
for (int ei = 0; ei < 3; ei++){
- int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3; //- iSign * iRotate
+ int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3;
+ //- iSign * iRotate
if (iSign > 0)
closure[fi + ei * (orderint - 1) + l] =
ti + edgenumber * (orderint - 1) + l;
@@ -832,8 +834,8 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
}
}
-static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<int> &closure,
- int order)
+static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
+ std::vector<int> &closure, int order)
{
if (order > 2)
Msg::Error("FaceClosure not implemented for prisms of order %d",order);
@@ -845,9 +847,12 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<i
closure[0] = 0;
return;
}
- int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
- int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15}, {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
-// int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8}, {8, 13, 11, 7}, {9, 11, 14, 10}};
+ int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2},
+ {1, 2, 5, 4}};
+ int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15},
+ {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
+ // int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8},
+ // {8, 13, 11, 7}, {9, 11, 14, 10}};
int nVertex = isTriangle ? 3 : 4;
for (int i = 0; i < nVertex; ++i){
int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
@@ -855,7 +860,8 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<i
}
if (order==2) {
for (int i = 0; i < nVertex; ++i){
- int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex;
+ //- iSign * iRotate
closure[nVertex+i] = order2node[iFace][k];
}
if (!isTriangle)
@@ -878,8 +884,8 @@ static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int ord
}
}
-
-static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
+static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order,
+ int nNod = 3)
{
closure.clear();
closure.resize(2*nNod);
@@ -897,13 +903,12 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, i
static void generate1dVertexClosure(polynomialBasis::clCont &closure)
{
-
closure.clear();
closure.resize(2);
closure[0].push_back(0);
closure[1].push_back(1);
-
}
+
std::map<int, polynomialBasis> polynomialBases::fs;
const polynomialBasis *polynomialBases::find(int tag)
@@ -1314,7 +1319,8 @@ std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
{
std::pair<int,int> key(tag1,tag2);
- std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
+ std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it =
+ injector.find(key);
if (it != injector.end()) return it->second;
const polynomialBasis& fs1 = *find(tag1);
@@ -1337,7 +1343,9 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
void polynomialBasis::registerBindings(binding *b) {
classBinding *cb = b->addClass<polynomialBasis>("polynomialBasis");
cb->setDescription("polynomial shape functions for elements");
- methodBinding *mb = cb->addMethod("f",(void (polynomialBasis::*)(fullMatrix<double>&, fullMatrix<double>&))&polynomialBasis::f);
+ methodBinding *mb = cb->addMethod
+ ("f", (void (polynomialBasis::*)(fullMatrix<double>&, fullMatrix<double>&))
+ &polynomialBasis::f);
mb->setDescription("evaluate the shape functions");
mb->setArgNames("nodes","values",NULL);
mb = cb->addMethod("find",&polynomialBases::find);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index c8038e2cc6..805df454d6 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -12,7 +12,8 @@
#include "fullMatrix.h"
#include <iostream>
-inline double pow_int (const double &a , const int &n) {
+inline double pow_int(const double &a, const int &n)
+{
switch (n) {
case 0 : return 1.0;
case 1 : return a;
@@ -77,11 +78,12 @@ class polynomialBasis
int numFaces;
// for a given face/edge, with both a sign and a rotation,
// give an ordered list of nodes on this face/edge
- inline const std::vector<int> &getClosure(int id) const // return the closure of dimension dim
+ inline const std::vector<int> &getClosure(int id) const
{
return closures[id];
}
- inline int getClosureId(int iEl, int iSign=1, int iRot=0) const {
+ inline int getClosureId(int iEl, int iSign=1, int iRot=0) const
+ {
return iEl + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
}
inline void evaluateMonomials(double u, double v, double w, double p[]) const
@@ -103,8 +105,10 @@ class polynomialBasis
}
}
}
- // I would favour this interface that allows optimizations (not implemented) and is easier to bind
- inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) {
+ // I would favour this interface that allows optimizations (not
+ // implemented) and is easier to bind
+ inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf)
+ {
double p[256];
sf.resize (coord.size1(), coefficients.size1());
for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index aba7c200c3..aca43d36f0 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -6,6 +6,7 @@
#ifndef _SIMPLE_FUNCTION_H_
#define _SIMPLE_FUNCTION_H_
+// FIXME: Numeric/ should not depend on Geo/
#include "MElement.h"
template <class scalar>
@@ -28,11 +29,12 @@ class simpleFunction {
#ifdef HAVE_LUA
#include "LuaBindings.h"
template <class scalar>
-class simpleFunctionLua:public simpleFunction<scalar> {
+class simpleFunctionLua : public simpleFunction<scalar> {
lua_State *_L;
std::string _luaFunctionName;
- public:
- scalar operator () (double x, double y, double z) const {
+ public:
+ scalar operator () (double x, double y, double z) const
+ {
lua_getfield(_L, LUA_GLOBALSINDEX, _luaFunctionName.c_str());
luaStack<double>::push(_L, x);
luaStack<double>::push(_L, y);
@@ -40,7 +42,9 @@ class simpleFunctionLua:public simpleFunction<scalar> {
lua_call(_L, 3, 1);
return luaStack<scalar>::get(_L,-1);
}
- simpleFunctionLua (lua_State *L, const std::string luaFunctionName, scalar s):simpleFunction<scalar>(s) {
+ simpleFunctionLua (lua_State *L, const std::string luaFunctionName, scalar s)
+ : simpleFunction<scalar>(s)
+ {
_L = L;
_luaFunctionName = luaFunctionName;
}
@@ -54,8 +58,8 @@ class simpleFunctionOnElement : public simpleFunction<scalar>
public :
simpleFunctionOnElement(scalar val=0) : simpleFunction<scalar>(val),_e(0) {}
virtual ~simpleFunctionOnElement(){}
- void setElement(MElement *e) {_e=e;}
- MElement * getElement(void) const {return _e;}
+ void setElement(MElement *e) { _e = e; }
+ MElement * getElement(void) const { return _e; }
MElement * getElement(double x, double y, double z) const
{
if (_e) return _e;
@@ -65,6 +69,4 @@ class simpleFunctionOnElement : public simpleFunction<scalar>
}
};
-
-
#endif
--
GitLab