From 8f230d42e3fb8442b346a8618ebdb34b7ae5a12d Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 21 Jul 2010 07:01:34 +0000
Subject: [PATCH] multilevel cartesian mostly done
---
Numeric/Numeric.cpp | 81 +++++----
Numeric/cartesian.h | 229 ++++++++++--------------
utils/api_demos/mainCartesian.cpp | 288 +++++++++++++++++++++---------
3 files changed, 346 insertions(+), 252 deletions(-)
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 73d6179585..b83bdf6ed3 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -754,39 +754,38 @@ double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
fullVector<double> grad(N);
fullVector<double> dir(N);
- double f,feps,finit;
+ double f, feps, finit;
for (int iter = 0; iter < MAXIT; iter++){
// compute gradient of func
- f = func(x,data);
- if (iter == 0)finit = f;
- // printf("Opti iter %d x = (%g %g) f = %g\n",iter,x(0),x(1),f);
- // printf("grad = (");
+ f = func(x, data);
+ if (iter == 0) finit = f;
+ // printf("Opti iter %d x = (%g %g) f = %g\n",iter,x(0),x(1),f);
+ // printf("grad = (");
for (int j = 0; j < N; j++){
double h = EPS * fabs(x(j));
if(h == 0.) h = EPS;
x(j) += h;
feps = func(x, data);
grad(j) = (feps - f) / h;
- // printf("%g ",grad(j));
+ // printf("%g ",grad(j));
dir(j) = -grad(j);
x(j) -= h;
}
- // printf(")\n ");
+ // printf(")\n ");
// do a 1D line search to fine the minimum
// of f(x - \alpha \nabla f)
double f, stpmax=100000;
int check;
- gmshLineSearch (func, data, x,dir,grad,f,stpmax,check);
- // printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
- if (check == 1)break;
+ gmshLineSearch(func, data, x, dir, grad, f, stpmax, check);
+ // printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
+ if (check == 1) break;
}
return f;
}
/*
-
P(p) = p1 + t1 xi + t2 eta
t1 = (p2-p1) ; t2 = (p3-p1) ;
@@ -809,12 +808,10 @@ distance to segment
- t ||p2-p1||^2 + (p-p1)(p2-p1) = 0
t = (p-p1)*(p2-p1)/||p2-p1||^2
-
-
*/
-void signedDistancesPointsTriangle(std::vector<double>&distances,
- std::vector<SPoint3>&closePts,
+void signedDistancesPointsTriangle(std::vector<double> &distances,
+ std::vector<SPoint3> &closePts,
const std::vector<SPoint3> &pts,
const SPoint3 &p1,
const SPoint3 &p2,
@@ -831,7 +828,7 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
{t1.z(), t2.z(), -n.z()}};
double inv[3][3];
double det = inv3x3(mat, inv);
- const unsigned pts_size=pts.size();
+ const unsigned pts_size = pts.size();
distances.clear();
distances.resize(pts_size);
closePts.clear();
@@ -847,7 +844,7 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
const double n2t3 = dot(t3, t3);
double u, v, d;
- for (unsigned int i = 0; i < pts_size;++i){
+ for (unsigned int i = 0; i < pts_size; ++i){
const SPoint3 &p = pts[i];
SVector3 pp1 = p - p1;
u = (inv[0][0] * pp1.x() + inv[0][1] * pp1.y() + inv[0][2] * pp1.z());
@@ -873,12 +870,12 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
found = true;
}
if (t13 >= 0 && t13 <= 1.){
- if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
+ if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
d = sign * std::min(fabs(d), p.distance(p1 + (p3 - p1) * t13));
found = true;
}
if (t23 >= 0 && t23 <= 1.){
- if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
+ if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
d = sign * std::min(fabs(d), p.distance(p2 + (p3 - p2) * t23));
found = true;
}
@@ -901,7 +898,9 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
}
}
-void signedDistancePointLine(const SPoint3 &p1,const SPoint3 &p2,const SPoint3 &p, double &d, SPoint3 &closePt){
+void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p,
+ double &d, SPoint3 &closePt)
+{
SVector3 t1 = p2 - p1;
const double n2t1 = dot(t1, t1);
SVector3 pp1 = p - p1;
@@ -923,12 +922,12 @@ void signedDistancePointLine(const SPoint3 &p1,const SPoint3 &p2,const SPoint3 &
}
}
-void signedDistancesPointsLine (std::vector<double>&distances,
- std::vector<SPoint3>&closePts,
- const std::vector<SPoint3> &pts,
- const SPoint3 &p1,
- const SPoint3 &p2){
-
+void signedDistancesPointsLine(std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2)
+{
distances.clear();
distances.resize(pts.size());
closePts.clear();
@@ -944,8 +943,11 @@ void signedDistancesPointsLine (std::vector<double>&distances,
}
}
-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){
- if (direction==1){
+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)
+{
+ if (direction == 1){
const SPoint3 &d1=SPoint3(1.0,0.0,0.0);
const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
@@ -962,7 +964,8 @@ void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &close
*x=closePt.x()*d1.x()+closePt.y()*d1.y()+closePt.z()*d1.z();
*y=closePt.x()*d3n.x()+closePt.y()*d3n.y()+closePt.z()*d3n.z();
*other=closePt.x()*d2n.x()+closePt.y()*d2n.y()+closePt.z()*d2n.z();
- }else{
+ }
+ else{
const SPoint3 &d2=SPoint3(0.0,1.0,0.0);
const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
@@ -982,7 +985,9 @@ void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &close
}
}
-int computeDistanceRatio(const double &y, const double &yp,const double &x,const double &xp, double *distance, const double &r1, const double &r2){
+int computeDistanceRatio(const double &y, const double &yp, const double &x,
+ const double &xp, double *distance, const double &r1, const double &r2)
+{
double b;
double a;
if (y==yp){
@@ -1083,14 +1088,14 @@ int computeDistanceRatio(const double &y, const double &yp,const double &x,const
}
}
-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){
-
+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)
+{
distances.clear();
distances.resize(pts.size());
distancesE.clear();
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index f064d24944..bd6f0d238e 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -11,111 +11,93 @@
// This is a cartesian mesh that encompasses an oriented box with NXI
// * NETA * NZETA hexaderal cells, with values stored at vertices
-/* j
+/*
+ 2-D example:
+
+ j
+---+---+---+---+---+---+
| | | |ij | | |
i +---+---+---+---+---+---+
| | | | | | |
+---+---+---+---+---+---+
- cell ij -> nodes i;j , i+i;j , i+1;j+1, i;j+1
-
- active = CELLS : store i+N*j
+ Nodal values/active cells are stored and referenced by node/cell
+ index, i.e., i+N*j
+
+ The ij cell has nodes i;j , i+1;j , i+1;j+1, i;j+1
*/
template <class scalar>
class cartesianBox {
private:
int _Nxi, _Neta, _Nzeta;
- std::set<int> _active;
double _X, _Y, _Z, _dxi, _deta, _dzeta;
SVector3 _xiAxis, _etaAxis, _zetaAxis;
+ std::set<int> _activeCells;
typename std::map<int, scalar> _nodalValues;
+ int _level;
+ cartesianBox<scalar> *_childBox;
public:
cartesianBox(double X, double Y, double Z,
const SVector3 &DXI, const SVector3 &DETA, const SVector3 &DZETA,
- int NXI, int NETA, int NZETA)
+ int NXI, int NETA, int NZETA, int level=1)
: _X(X), _Y(Y), _Z(Z),
_dxi(norm(DXI)), _deta(norm(DETA)), _dzeta(norm(DZETA)),
_xiAxis(DXI), _etaAxis(DETA), _zetaAxis(DZETA),
- _Nxi(NXI), _Neta(NETA), _Nzeta(NZETA)
+ _Nxi(NXI), _Neta(NETA), _Nzeta(NZETA), _level(level), _childBox(0)
{
_xiAxis.normalize();
_etaAxis.normalize();
_zetaAxis.normalize();
+ if(level > 1)
+ _childBox = new cartesianBox<scalar>(X, Y, Z, DXI, DETA, DZETA,
+ 2 * NXI, 2 * NETA, 2 * NZETA,
+ level - 1);
}
- typename std::map<int, scalar>::const_iterator begin() const
- {
- return _nodalValues.begin();
- }
- typename std::map<int, scalar>::const_iterator end() const
- {
- return _nodalValues.end();
- }
- typedef std::set<int>::iterator boxIter;
- boxIter activeBoxBegin() { return _active.begin(); }
- boxIter activeBoxEnd() { return _active.end(); }
+ int getNxi(){ return _Nxi; }
+ int getNeta(){ return _Neta; }
+ int getNzeta(){ return _Nzeta; }
+ cartesianBox<scalar> *getChildBox(){ return _childBox; }
+ int getLevel(){ return _level; }
+ typedef std::set<int>::const_iterator cellIter;
+ cellIter activeCellsBegin(){ return _activeCells.begin(); }
+ cellIter activeCellsEnd(){ return _activeCells.end(); }
+ typedef typename std::map<int, scalar>::iterator valIter;
+ valIter nodalValuesBegin(){ return _nodalValues.begin(); }
+ valIter nodalValuesEnd(){ return _nodalValues.end(); }
+ void setNodalValue(int i, scalar s){ _nodalValues[i] = s; }
void getNodalValues(const int &t, std::vector<scalar> &ls_values)
{
// perhaps not optimal (lot of searches...)
int i, j, k;
- element_ijk(t, i, j, k);
+ getCellIJK(t, i, j, k);
typename std::map<int, scalar>::iterator itNode;
- itNode = _nodalValues.find(node_index(i, j, k));
+ itNode = _nodalValues.find(getNodeIndex(i, j, k));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i+1,j,k));
+ itNode = _nodalValues.find(getNodeIndex(i+1,j,k));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i+1,j+1,k));
+ itNode = _nodalValues.find(getNodeIndex(i+1,j+1,k));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i,j+1,k));
+ itNode = _nodalValues.find(getNodeIndex(i,j+1,k));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i,j,k+1));
+ itNode = _nodalValues.find(getNodeIndex(i,j,k+1));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i+1,j,k+1));
+ itNode = _nodalValues.find(getNodeIndex(i+1,j,k+1));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i+1,j+1,k+1));
+ itNode = _nodalValues.find(getNodeIndex(i+1,j+1,k+1));
ls_values.push_back(itNode->second);
- itNode = _nodalValues.find(node_index(i,j+1,k+1));
+ itNode = _nodalValues.find(getNodeIndex(i,j+1,k+1));
ls_values.push_back(itNode->second);
}
-
-// // Test : does not work for the moment
-// bool connectedToActiveCell(int i, int j, int k){
-//
-// int connectedCells[8]={-1,-1,-1,-1,-1,-1,-1,-1};
-//
-// connectedCells[0]=element_index(i,j,k);
-// connectedCells[1]=element_index(i-1,j,k);
-// connectedCells[2]=element_index(i,j-1,k);
-// connectedCells[3]=element_index(i-1,j-1,k);
-// connectedCells[4]=element_index(i,j,k-1);
-// connectedCells[5]=element_index(i-1,j,k-1);
-// connectedCells[6]=element_index(i,j-1,k-1);
-// connectedCells[7]=element_index(i-1,j-1,k-1);
-//
-// bool activeCell=false;
-//
-// for(int I=0;I<8;++I){
-// activeCell += (_active.find() == _active.end());
-// }
-//
-// return activeCell;
-//
-// }
- // set the value
- void setValue(int i, scalar s)
- {
- _nodalValues[i] = s;
- }
- inline int index_of_element(double x, double y, double z) const
+ inline int getCellContainingPoint(double x, double y, double z) const
{
// P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
// DP = P-P_0 * _vdx = xi
@@ -123,79 +105,64 @@ class cartesianBox {
double xi = dot(DP, _xiAxis);
double eta = dot(DP, _etaAxis);
double zeta = dot(DP, _zetaAxis);
-
int i = xi / _dxi * _Nxi;
int j = eta / _deta * _Neta;
int k = zeta / _dzeta * _Nzeta;
-
- if (i < 0) i = 0;if (i >= _Nxi) i = _Nxi - 1;
- if (j < 0) j = 0;if (j >= _Neta) j = _Neta - 1;
- if (k < 0) k = 0;if (k >= _Nzeta) k = _Nzeta - 1;
-
- return element_index(i, j, k);
+ if (i < 0) i = 0; if (i >= _Nxi) i = _Nxi - 1;
+ if (j < 0) j = 0; if (j >= _Neta) j = _Neta - 1;
+ if (k < 0) k = 0; if (k >= _Nzeta) k = _Nzeta - 1;
+ return getCellIndex(i, j, k);
}
- inline SPoint3 coordinates_of_node(const int &t) const
+ inline SPoint3 getNodeCoordinates(const int &t) const
{
int i, j, k;
- node_ijk(t, i, j, k);
-
- // SVector3 DP (x-_X,y-_Y,z-_Z);
-
+ getNodeIJK(t, i, j, k);
const double xi = i * _dxi / _Nxi;
const double eta = j * _deta / _Neta;
const double zeta = k * _dzeta / _Nzeta;
-
- // printf("index %d -> %d %d %d xi %g %g %g X %g %g %g N %d %d %d D %g %g %g\n",
- // t,i,j,k,xi,eta,zeta,_X,_Y,_Z,_Nxi,_Neta,_Nzeta,_dxi,_deta,_dzeta);
-
SVector3 D = xi * _xiAxis + eta * _etaAxis + zeta * _zetaAxis;
-
return SPoint3(_X + D.x(), _Y + D.y(), _Z + D.z());
}
- void insert(const int &t)
- {
- _active.insert(t);
- }
- void erase(const int &t)
- {
- _active.erase(t);
- }
- void erase(std::set<int>::iterator t)
+ void insertActiveCell(const int &t){ _activeCells.insert(t); }
+ void eraseActiveCell(const int &t){ _activeCells.erase(t); }
+ void eraseActiveCell(std::set<int>::iterator t){ _activeCells.erase(t); }
+ bool activeCellExists(const int &t)
{
- _active.erase(t);
+ return (_activeCells.find(t) != _activeCells.end());
}
- inline int element_index(int i, int j, int k) const
+ inline int getCellIndex(int i, int j, int k) const
{
return i + _Nxi * j + _Nxi *_Neta * k;
}
- inline int node_index(int i, int j, int k) const
+ inline int getNodeIndex(int i, int j, int k) const
{
return i + (_Nxi+1) * j + (_Nxi+1) * (_Neta+1) * k;
}
- inline void element_ijk(int index, int &i, int &j, int &k) const
+ inline void getCellIJK(int index, int &i, int &j, int &k) const
{
k = index / (_Nxi * _Neta);
j = (index - k * (_Nxi * _Neta)) / _Nxi;
i = (index - k * (_Nxi * _Neta) - j * _Nxi);
}
- inline void node_ijk(int index, int &i, int &j, int &k) const
+ inline void getNodeIJK(int index, int &i, int &j, int &k) const
{
k = index / ((_Nxi + 1) * (_Neta + 1));
j = (index - k * ((_Nxi + 1) * (_Neta + 1))) / (_Nxi + 1);
i = (index - k * ((_Nxi + 1) * (_Neta + 1)) - j * (_Nxi + 1));
}
- inline void create_nodes()
+ inline void createNodalValues()
{
- std::set<int>::const_iterator it = _active.begin();
- for( ; it != _active.end() ; ++it){
+ std::set<int>::const_iterator it = _activeCells.begin();
+ for( ; it != _activeCells.end() ; ++it){
const int &t = *it;
int i, j, k;
- element_ijk(t, i, j, k);
+ getCellIJK(t, i, j, k);
for (int I = 0; I < 2; I++)
for (int J = 0; J < 2; J++)
for (int K = 0; K < 2; K++)
- _nodalValues[node_index(i + I, j + J, k + K)] = 0.;
+ _nodalValues[getNodeIndex(i + I, j + J, k + K)] = 0.;
}
+ if(_childBox) _childBox->createNodalValues();
}
void writeMSH(const std::string &filename, bool simplex=false,
bool writeNodalValues=true) const
@@ -206,76 +173,78 @@ class cartesianBox {
fprintf(f, "$Nodes\n%d\n", (int)_nodalValues.size());
typename std::map<int, scalar>::const_iterator it = _nodalValues.begin();
for ( ; it != _nodalValues.end(); ++it){
- SPoint3 p = coordinates_of_node(it->first);
+ SPoint3 p = getNodeCoordinates(it->first);
fprintf(f, "%d %g %g %g\n", it->first, p.x(), p.y(), p.z());
}
fprintf(f, "$EndNodes\n");
}
{
int coeff = simplex ? 6 : 1;
- fprintf(f,"$Elements\n%d\n", coeff * (int)_active.size());
- std::set<int>::const_iterator it = _active.begin();
- for ( ; it != _active.end(); ++it){
+ fprintf(f,"$Elements\n%d\n", coeff * (int)_activeCells.size());
+ std::set<int>::const_iterator it = _activeCells.begin();
+ for ( ; it != _activeCells.end(); ++it){
int i, j, k;
- element_ijk(*it, i, j, k);
+ getCellIJK(*it, i, j, k);
if(!simplex){
fprintf(f, "%d 5 3 1 1 1", *it);
- fprintf(f, " %d", node_index(i, j, k));
- fprintf(f, " %d", node_index(i + 1, j, k));
- fprintf(f, " %d", node_index(i + 1, j + 1, k));
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j + 1, k + 1));
- fprintf(f, " %d", node_index(i, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k + 1));
fprintf(f, "\n");
}
else{
// Elt1
fprintf(f, "%d 4 3 1 1 1", *it);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i, j + 1, k + 1));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k + 1));
fprintf(f, "\n");
// Elt2
fprintf(f, "%d 4 3 1 1 1", *it + 1);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i + 1, j + 1, k + 1));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k));
fprintf(f, "\n");
// Elt3
fprintf(f, "%d 4 3 1 1 1", *it + 2);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i, j + 1, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k + 1));
fprintf(f, "\n");
// Elt4
fprintf(f, "%d 4 3 1 1 1", *it + 3);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i + 1, j + 1, k));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i + 1, j, k));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k));
fprintf(f, "\n");
// Elt5
fprintf(f, "%d 4 3 1 1 1", *it + 4);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i + 1, j, k));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i, j, k));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j, k));
fprintf(f, "\n");
// Elt6
fprintf(f, "%d 4 3 1 1 1", *it + 5);
- fprintf(f, " %d", node_index(i, j + 1, k));
- fprintf(f, " %d", node_index(i, j, k));
- fprintf(f, " %d", node_index(i + 1, j, k + 1));
- fprintf(f, " %d", node_index(i, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j + 1, k));
+ fprintf(f, " %d", getNodeIndex(i, j, k));
+ fprintf(f, " %d", getNodeIndex(i + 1, j, k + 1));
+ fprintf(f, " %d", getNodeIndex(i, j, k + 1));
fprintf(f, "\n");
}
}
fprintf(f,"$EndElements\n");
+ if(_childBox)
+ _childBox->writeMSH(filename + "_2", simplex, writeNodalValues);
}
if(writeNodalValues){
fprintf(f,"$NodeData\n1\n\"distance\"\n 1\n 0.0\n3\n0\n 1\n %d\n",
diff --git a/utils/api_demos/mainCartesian.cpp b/utils/api_demos/mainCartesian.cpp
index de3f52327d..052ce22e35 100644
--- a/utils/api_demos/mainCartesian.cpp
+++ b/utils/api_demos/mainCartesian.cpp
@@ -5,125 +5,246 @@
#include "MTriangle.h"
#include "SOrientedBoundingBox.h"
#include "Numeric.h"
+#include "GmshMessage.h"
+#include "OS.h"
-void insertBoxes(double x, double y, double z, double thickness, cartesianBox<double> &box)
+void insertActiveCells(double x, double y, double z, double rmax,
+ cartesianBox<double> &box)
{
- int id1 = box.index_of_element(x - thickness, y - thickness, z - thickness);
- int id2 = box.index_of_element(x + thickness, y + thickness, z + thickness);
+ int id1 = box.getCellContainingPoint(x - rmax, y - rmax, z - rmax);
+ int id2 = box.getCellContainingPoint(x + rmax, y + rmax, z + rmax);
int i1, j1 ,k1;
- box.element_ijk(id1, i1, j1, k1);
+ box.getCellIJK(id1, i1, j1, k1);
int i2, j2, k2;
- box.element_ijk(id2, i2, j2, k2);
- for (int i = i1; i <= i2; i++){
- for (int j = j1; j <= j2; j++){
- for (int k = k1; k <= k2; k++){
- int id = box.element_index(i, j, k);
- box.insert(id);
- }
+ box.getCellIJK(id2, i2, j2, k2);
+ for (int i = i1; i <= i2; i++)
+ for (int j = j1; j <= j2; j++)
+ for (int k = k1; k <= k2; k++)
+ box.insertActiveCell(box.getCellIndex(i, j, k));
+}
+
+void computeLevelset(GModel *gm, cartesianBox<double> &box)
+{
+ std::vector<SPoint3> nodes;
+ std::vector<int> indices;
+ for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
+ it != box.nodalValuesEnd(); ++it){
+ nodes.push_back(box.getNodeCoordinates(it->first));
+ indices.push_back(it->first);
+ }
+ Msg::Info(" %d nodes in the grid", (int)nodes.size());
+ std::vector<double> dist, localdist;
+ std::vector<SPoint3> dummy;
+ for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
+ for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
+ int i1 = (*fit)->stl_triangles[i];
+ int i2 = (*fit)->stl_triangles[i + 1];
+ int i3 = (*fit)->stl_triangles[i + 2];
+ GPoint p1 = (*fit)->point((*fit)->stl_vertices[i1]);
+ GPoint p2 = (*fit)->point((*fit)->stl_vertices[i2]);
+ GPoint p3 = (*fit)->point((*fit)->stl_vertices[i3]);
+ SPoint2 p = ((*fit)->stl_vertices[i1] + (*fit)->stl_vertices[i2] +
+ (*fit)->stl_vertices[i3]) * 0.33333333;
+ SVector3 N = (*fit)->normal(p);
+ SPoint3 P1(p1.x(), p1.y(), p1.z());
+ SPoint3 P2(p2.x(), p2.y(), p2.z());
+ SPoint3 P3(p3.x(), p3.y(), p3.z());
+ SVector3 NN(crossprod(P2 - P1, P3 - P1));
+ if (dot(NN, N) > 0)
+ signedDistancesPointsTriangle(localdist, dummy, nodes, P1, P2, P3);
+ else
+ signedDistancesPointsTriangle(localdist, dummy, nodes, P2, P1, P3);
+ if(dist.empty())
+ dist = localdist;
+ else
+ for (unsigned int j = 0; j < localdist.size(); j++)
+ dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
}
}
+ for (unsigned int j = 0; j < dist.size(); j++)
+ box.setNodalValue(indices[j], dist[j]);
+
+ if(box.getChildBox()) computeLevelset(gm, *box.getChildBox());
+}
+
+void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
+{
+ Range<double> t_bounds = ge->parBounds(0);
+ double t_min = t_bounds.low();
+ double t_max = t_bounds.high();
+ double length = ge->length(t_min, t_max, 20);
+ int N = length / sampling;
+ for(int i = 0; i < N; i++) {
+ double t = t_min + (double)i / (double)(N - 1) * (t_max - t_min);
+ GPoint p = ge->point(t);
+ double x = p.x(), y = p.y(), z = p.z();
+ points.push_back(SPoint3(x, y, z));
+ }
+}
+
+int removeOrphanChildCells(cartesianBox<double> *parent)
+{
+ cartesianBox<double> *child = parent->getChildBox();
+ if(!child) return 0;
+ int I = parent->getNxi();
+ int J = parent->getNeta();
+ int K = parent->getNzeta();
+ for(int i = 0; i < I; i++)
+ for(int j = 0; j < J; j++)
+ for(int k = 0; k < K; k++){
+ // remove children if they do not entirely fill parent, or if
+ // there is no parent on the boundary
+ int idx[8] = {child->getCellIndex(2 * i, 2 * j, 2 * k),
+ child->getCellIndex(2 * i, 2 * j, 2 * k + 1),
+ child->getCellIndex(2 * i, 2 * j + 1, 2 * k),
+ child->getCellIndex(2 * i, 2 * j + 1, 2 * k + 1),
+ child->getCellIndex(2 * i + 1, 2 * j, 2 * k),
+ child->getCellIndex(2 * i + 1, 2 * j, 2 * k + 1),
+ child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k),
+ child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k + 1)};
+ bool exist[8], atLeastOne = false, butNotAll = false;
+ for(int ii = 0; ii < 8; ii++){
+ exist[ii] = child->activeCellExists(idx[ii]);
+ if(exist[ii]) atLeastOne = true;
+ if(!exist[ii]) butNotAll = true;
+ }
+ if(atLeastOne &&
+ (butNotAll ||
+ (i != 0 && !parent->activeCellExists(parent->getCellIndex(i - 1, j, k))) ||
+ (i != I - 1 && !parent->activeCellExists(parent->getCellIndex(i + 1, j, k))) ||
+ (j != 0 && !parent->activeCellExists(parent->getCellIndex(i, j - 1, k))) ||
+ (j != J - 1 && !parent->activeCellExists(parent->getCellIndex(i, j + 1, k))) ||
+ (k != 0 && !parent->activeCellExists(parent->getCellIndex(i, j, k - 1))) ||
+ (k != K - 1 && !parent->activeCellExists(parent->getCellIndex(i, j, k + 1)))))
+ for(int ii = 0; ii < 8; ii++) child->eraseActiveCell(idx[ii]);
+ }
+ return removeOrphanChildCells(child);
}
-int main (int argc,char *argv[])
+void removeParentCellsWithChildren(cartesianBox<double> *box)
{
- if(argc < 4){
- printf("Usage: %s file thickness npointsx [sampling] [filter]\n", argv[0]);
+ if(!box->getChildBox()) return;
+ for(int i = 0; i < box->getNxi(); i++)
+ for(int j = 0; j < box->getNeta(); j++)
+ for(int k = 0; k < box->getNzeta(); k++){
+ if(box->activeCellExists(box->getCellIndex(i, j, k))){
+ cartesianBox<double> *parent = box, *child;
+ int ii = i, jj = j, kk = k;
+ while((child = parent->getChildBox())){
+ ii *= 2; jj *= 2; kk *= 2;
+ if(child->activeCellExists(child->getCellIndex(ii, jj, kk))){
+ box->eraseActiveCell(box->getCellIndex(i, j, k));
+ break;
+ }
+ parent = child;
+ }
+ }
+ }
+ removeParentCellsWithChildren(box->getChildBox());
+}
+
+int main(int argc,char *argv[])
+{
+ if(argc < 6){
+ printf("Usage: %s file lx ly lz rmax [levels=1] [filter]\n", argv[0]);
printf("where\n");
- printf(" 'file' contains a CAD model or a mesh\n");
- printf(" 'thickness' sets the thickness of the layer of elements created around the surfaces\n");
- printf(" 'npointsx' sets the number of mesh points along the x-axis\n");
- printf(" 'sampling' sets the maximum distance between 2 sampling points on the surfaces\n");
- printf(" 'filter' selects if only cells inside the body are saved in the mesh\n");
+ printf(" 'file' contains a CAD model\n");
+ printf(" 'lx', 'ly' and 'lz' are the sizes of the elements along the"
+ " x-, y- and z-axis at the coarsest level\n");
+ printf(" 'rmax' is the radius of the largest sphere that can be inscribed"
+ " in the structure\n");
+ printf(" 'levels' sets the number of levels in the grid\n");
+ printf(" 'filter' selects if only cells inside the body are saved in the"
+ " mesh\n");
return -1;
}
GmshInitialize();
+ GmshSetOption("General", "Terminal", 1.);
GmshMergeFile(argv[1]);
- double thickness = atof(argv[2]);
- int npointsx = atoi(argv[3]);
- double sampling = (argc > 4) ? atof(argv[4]) : thickness / 3.;
- int filter = (argc > 5) ? atoi(argv[5]) : 0;
+ double lx = atof(argv[2]), ly = atof(argv[3]), lz = atof(argv[4]);
+ double rmax = atof(argv[5]);
+ int levels = (argc > 6) ? atof(argv[6]) : 1;
+ int filter = (argc > 7) ? atoi(argv[7]) : 0;
+
+ // minimum distance between points in the cloud at the coarsest level
+ double sampling = std::min(lx, std::min(ly, lz)) / 2.;
+
+ // radius of the "tube" created around edges at the coarsest level
+ double rtube = std::max(lx, std::max(ly, lz)) * 2.;
GModel *gm = GModel::current();
std::vector<SPoint3> points;
- std::vector<SVector3> normals;
-
+ Msg::Info("Filling point clould on surfaces");
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
- (*fit)->fillPointCloud(sampling, &points, &normals);
- printf("%d points in the cloud\n", (int)points.size());
+ (*fit)->fillPointCloud(sampling, &points);
+ Msg::Info(" %d points in the surface cloud", (int)points.size());
- // printf("OBB with %d points\n", points.size());
- // SOrientedBoundingBox *obb = SOrientedBoundingBox::buildOBB(points);
- // printf("OBB done\n");
+ std::vector<SPoint3> edgePoints;
+ if(levels > 1){
+ Msg::Info("Filling point clould on curves");
+ for (GModel::eiter eit = gm->firstEdge(); eit != gm->lastEdge(); eit++)
+ fillPointCloud(*eit, sampling / pow(2., levels - 1), edgePoints);
+ Msg::Info(" %d points in the curvee cloud", (int)edgePoints.size());
+ }
SBoundingBox3d bb;
for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
+ for(unsigned int i = 0; i < edgePoints.size(); i++) bb += edgePoints[i];
bb.scale(1.2, 1.2, 1.2);
SVector3 range = bb.max() - bb.min();
- int NX = npointsx;
- int NY = NX * range.y() / range.x();
- int NZ = NX * range.z() / range.x();
+ int NX = range.x() / lx;
+ int NY = range.y() / ly;
+ int NZ = range.z() / lz;
+ if(NX < 2) NX = 2;
+ if(NY < 2) NY = 2;
+ if(NZ < 2) NZ = 2;
- printf("bb Min= %g %g %g -- bb Max= %g %g %g --NX %d NY %d NZ %d\n",
- bb.min().x(), bb.min().y(), bb.min().z(),
- bb.max().x(), bb.max().y(), bb.max().z(), NX, NY, NZ);
+ Msg::Info(" bb Min= %g %g %g -- bb Max= %g %g %g -- NX %d NY %d NZ %d",
+ bb.min().x(), bb.min().y(), bb.min().z(),
+ bb.max().x(), bb.max().y(), bb.max().z(), NX, NY, NZ);
cartesianBox<double> box(bb.min().x(), bb.min().y(), bb.min().z(),
SVector3(range.x(), 0, 0),
SVector3(0, range.y(), 0),
SVector3(0, 0, range.z()),
- NX, NY, NZ);
+ NX, NY, NZ, levels);
+ Msg::Info("Inserting active cells in the cartesian grid");
for (unsigned int i = 0; i < points.size();i++)
- insertBoxes(points[i].x(), points[i].y(), points[i].z(), thickness, box);
+ insertActiveCells(points[i].x(), points[i].y(), points[i].z(), rmax, box);
+ cartesianBox<double> *parent = &box, *child;
+ while((child = parent->getChildBox())){
+ for(unsigned int i = 0; i < edgePoints.size(); i++)
+ insertActiveCells(edgePoints[i].x(), edgePoints[i].y(), edgePoints[i].z(),
+ rtube / pow(2., levels - child->getLevel()), *child);
+ parent = child;
+ }
- printf("insertion done\n");
- box.create_nodes();
+ Msg::Info("Removing cells to match X-FEM mesh topology constraints");
+ // remove child cells that do not entirely fill parent cell or for
+ // which there is no parent neighbor
+ removeOrphanChildCells(&box);
+ // remove parent cells that have children
+ removeParentCellsWithChildren(&box);
- std::map<int,double>::const_iterator it = box.begin();
- std::vector<SPoint3> NODES;
- std::vector<SPoint3> CNODES;
- std::vector<int> indices;
- std::vector<double> dist, localdist;
- for ( ; it != box.end(); ++it){
- NODES.push_back(box.coordinates_of_node(it->first));
- indices.push_back(it->first);
- }
- printf("%d nodes in the cartesian mesh\n", (int)NODES.size());
+ // TODO:
+ // * remove child nodes which exist in parent grid
+ // * offset of node numbers depending on box level to write a single mesh
- for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
- for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
- int i1 = (*fit)->stl_triangles[i];
- int i2 = (*fit)->stl_triangles[i + 1];
- int i3 = (*fit)->stl_triangles[i + 2];
- GPoint p1 = (*fit)->point((*fit)->stl_vertices[i1]);
- GPoint p2 = (*fit)->point((*fit)->stl_vertices[i2]);
- GPoint p3 = (*fit)->point((*fit)->stl_vertices[i3]);
- SPoint2 ppp = ((*fit)->stl_vertices[i1] + (*fit)->stl_vertices[i2] + (*fit)->stl_vertices[i3]) * 0.33333333;
- SVector3 N = (*fit)->normal(ppp);
- SPoint3 P1(p1.x(), p1.y(), p1.z());
- SPoint3 P2(p2.x(), p2.y(), p2.z());
- SPoint3 P3(p3.x(), p3.y(), p3.z());
- SVector3 NN(crossprod(P2 - P1, P3 - P1));
- if (dot(NN, N) > 0)
- signedDistancesPointsTriangle(localdist, CNODES, NODES, P1, P2, P3);
- else
- signedDistancesPointsTriangle(localdist, CNODES, NODES, P2, P1, P3);
- if(dist.empty())
- dist = localdist;
- else
- for (unsigned int j = 0; j < localdist.size(); j++)
- dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
- }
- }
- for (unsigned int j = 0; j < dist.size(); j++)
- box.setValue(indices[j], dist[j]);
+ Msg::Info("Initializing nodal values in the cartesian grid");
+ box.createNodalValues();
+
+ Msg::Info("Computing levelset on the cartesian grid");
+ computeLevelset(gm, box);
if(filter){
+ Msg::Info("Filtering result -- *** TODO THIS NEEDS TO BE ADAPTED FOR MULTILEVEL");
int nbErased = 0;
//Coup de menage avant d'exporter le maillage
- for(cartesianBox<double>::boxIter it = box.activeBoxBegin(); it != box.activeBoxEnd();){
+ for(cartesianBox<double>::cellIter it = box.activeCellsBegin();
+ it != box.activeCellsEnd();){
std::vector<double> ls_vals;
box.getNodalValues(*it, ls_vals);
double lsmax = *std::max_element(ls_vals.begin(), ls_vals.end());
@@ -131,19 +252,18 @@ int main (int argc,char *argv[])
double change_sign = lsmax * lsmin ;
double epsilon = 1.e-10;
if(change_sign > 0 && lsmax < -epsilon) {
- box.erase(*(it++));
+ box.eraseActiveCell(*(it++));
nbErased++;
}
else ++it;
}
- std::cout << "Number of erased cells after filtering : " << nbErased << std::endl;
+ Msg::Info(" number of erased cells after filtering: %d", nbErased);
}
- printf("nodes created\n");
- box.writeMSH("yeah.msh", true);
-
+ Msg::Info("Writing results to disk");
+ box.writeMSH("yeah.msh");
box.writeMSH("youhou.msh", true, false);
box.writeNodalValues("youhou.pos");
- printf("mesh written\n");
+ Msg::Info("Done!");
GmshFinalize();
}
--
GitLab