From 9df7372be05c41f09c6e467e9ebe07d0bd86b111 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 16 Jul 2010 13:14:11 +0000
Subject: [PATCH] cleanup cartesian code before more invasive changes for
multi-level grid
---
Numeric/CMakeLists.txt | 1 -
Numeric/cartesian.cpp | 52 -----
Numeric/cartesian.h | 327 +++++++++++++-----------------
utils/api_demos/mainCartesian.cpp | 258 ++++++++---------------
4 files changed, 229 insertions(+), 409 deletions(-)
delete mode 100644 Numeric/cartesian.cpp
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 4d7f5c9233..3a9b71fe38 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -19,7 +19,6 @@ set(SRC
mathEvaluator.cpp
Iso.cpp
DivideAndConquer.cpp
- cartesian.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Numeric/cartesian.cpp b/Numeric/cartesian.cpp
deleted file mode 100644
index 8b6193d51d..0000000000
--- a/Numeric/cartesian.cpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#include "GmshConfig.h"
-#include "cartesian.h"
-
-#if defined(HAVE_ANN)
-#include "ANN/ANN.h"
-class hiddenANN{
-public:
- ANNkd_tree *kdtree;
- ANNpointArray zeronodes;
- ANNidxArray index;
- ANNdistArray dist;
-};
-
-template <>
-double cartesianBox<double>::distance (double x, double y, double z) const
-{
- if (!_ann){
- _ann = new hiddenANN;
- _ann->index = new ANNidx[2];
- _ann->dist = new ANNdist[2];
- if(_points.size())
- _ann->zeronodes = annAllocPts(_points.size(), 4);
- for (unsigned int i = 0; i < _points.size(); i++){
- _ann->zeronodes[i][0] = _points[i].x();
- _ann->zeronodes[i][1] = _points[i].y();
- _ann->zeronodes[i][2] = _points[i].z();
- }
- _ann->kdtree = new ANNkd_tree(_ann->zeronodes, _points.size(), 3);
- }
- double xyz[3] = { x, y, z };
- _ann->kdtree->annkSearch(xyz, 1, _ann->index, _ann->dist);
- SVector3 n1 = _normals[_ann->index[0]];
- // SVector3 n2 = _normals[_ann->index[1]];
- SVector3 dx1 = SPoint3(x,y,z) - _points[_ann->index[0]];
- // SVector3 dx2 = SPoint3(x,y,z) - _points[_ann->index[1]];
- double d1 = sqrt(_ann->dist[0]);
- // double d2 = sqrt(_ann->dist[1]);
- double sign1 = dot(n1,dx1) > 0 ? 1. : -1.0;
- // double sign2 = dot(n2,dx2) > 0 ? 1. : -1.0;
- // if (fabs(d1-d2) > 1.e-4 * _dxi/((double)_Nxi)){
- return sign1 * d1;
- // }
- // return (sign1 != sign2) ? d1 : sign1 *d1;
-}
-#else
-template <>
-double cartesianBox<double>::distance (double x, double y, double z) const
-{
- printf("youpiiie\n");
- return 0.0;
-}
-#endif
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index a3285dd581..f064d24944 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -8,8 +8,8 @@
#include "SVector3.h"
#include "SPoint3.h"
-// this is a cartesian mesh that encompasses an oriented box with _NX
-// x _NY x _NZ hexaderal cells
+// This is a cartesian mesh that encompasses an oriented box with NXI
+// * NETA * NZETA hexaderal cells, with values stored at vertices
/* j
+---+---+---+---+---+---+
@@ -23,54 +23,46 @@
active = CELLS : store i+N*j
*/
-class hiddenANN;
-
template <class scalar>
class cartesianBox {
- mutable hiddenANN *_ann;
private:
int _Nxi, _Neta, _Nzeta;
std::set<int> _active;
double _X, _Y, _Z, _dxi, _deta, _dzeta;
SVector3 _xiAxis, _etaAxis, _zetaAxis;
typename std::map<int, scalar> _nodalValues;
- std::vector<SVector3> _normals;
- std::vector<SPoint3> _points;
public:
- std::vector<SPoint3> & points() { return _points; }
- std::vector<SVector3> & normals() { return _normals; }
- cartesianBox (double X, double Y, double Z,
- const SVector3 &DXI,
- const SVector3 &DETA,
- const SVector3 &DZETA,
- int NXI, int NETA, int NZETA)
- : _ann(0), _X(X), _Y(Y), _Z(Z),
- _dxi(norm(DXI)),
- _deta(norm(DETA)),
- _dzeta(norm(DZETA)),
- _xiAxis(DXI),
- _etaAxis(DETA),
- _zetaAxis(DZETA),
+ cartesianBox(double X, double Y, double Z,
+ const SVector3 &DXI, const SVector3 &DETA, const SVector3 &DZETA,
+ int NXI, int NETA, int NZETA)
+ : _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)
{
_xiAxis.normalize();
_etaAxis.normalize();
_zetaAxis.normalize();
}
-
- typename std::map<int, scalar>::const_iterator begin() const { return _nodalValues.begin(); }
- typename std::map<int, scalar>::const_iterator end() const { return _nodalValues.end(); }
-
+ 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(); }
-
- //Perhaps not optimal (lot of searches...)
- void getNodalValues(const int &t, std::vector<scalar> &ls_values) {
- int i,j,k;
- element_ijk(t,i,j,k);
+ boxIter activeBoxBegin() { return _active.begin(); }
+ boxIter activeBoxEnd() { return _active.end(); }
+ 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);
+ typename std::map<int, scalar>::iterator itNode;
- typename std::map<int, scalar>::iterator itNode = _nodalValues.find(node_index(i,j,k));
+ itNode = _nodalValues.find(node_index(i, j, k));
ls_values.push_back(itNode->second);
itNode = _nodalValues.find(node_index(i+1,j,k));
@@ -93,8 +85,6 @@ class cartesianBox {
itNode = _nodalValues.find(node_index(i,j+1,k+1));
ls_values.push_back(itNode->second);
-
- return;
}
// // Test : does not work for the moment
@@ -120,225 +110,192 @@ class cartesianBox {
// return activeCell;
//
// }
-
-
- // add that in the ann search tool
- void insert_point (double x, double y, double z)
- {
- _points.push_back(SPoint3(x,y,z));
- }
- // compute distance
- double distance (double x, double y, double z) const;
// 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 index_of_element(double x, double y, double z) const
{
// P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
// DP = P-P_0 * _vdx = xi
- SVector3 DP (x-_X,y-_Y,z-_Z);
+ SVector3 DP (x - _X, y - _Y, z - _Z);
+ double xi = dot(DP, _xiAxis);
+ double eta = dot(DP, _etaAxis);
+ double zeta = dot(DP, _zetaAxis);
- 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;
- 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;
- 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);
+ return element_index(i, j, k);
}
- inline SPoint3 coordinates_of_node (const int &t) const
+ inline SPoint3 coordinates_of_node(const int &t) const
{
int i, j, k;
node_ijk(t, i, j, k);
// SVector3 DP (x-_X,y-_Y,z-_Z);
- const double xi = i*_dxi/_Nxi;
- const double eta = j*_deta/_Neta;
- const double zeta = k*_dzeta/_Nzeta;
+ 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;
+ SVector3 D = xi * _xiAxis + eta * _etaAxis + zeta * _zetaAxis;
- return SPoint3(_X+D.x(),_Y+D.y(),_Z+D.z());
+ return SPoint3(_X + D.x(), _Y + D.y(), _Z + D.z());
}
- void insert (const int &t)
+ void insert(const int &t)
{
_active.insert(t);
}
- void erase (const int &t)
+ void erase(const int &t)
{
_active.erase(t);
}
-
- void erase (std::set<int>::iterator t)
+ void erase(std::set<int>::iterator t)
{
_active.erase(t);
}
-
- inline int element_index (int i, int j, int k) const
+ inline int element_index(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 node_index(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 element_ijk(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);
+ 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 node_ijk(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));
+ 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()
{
std::set<int>::const_iterator it = _active.begin();
for( ; it != _active.end() ; ++it){
const int &t = *it;
- int i,j,k;
- element_ijk(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.0;
+ int i, j, k;
+ element_ijk(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.;
}
}
- void writeMSH (const std::string &filename, bool simplex=false, bool writeLS=true) const
+ void writeMSH(const std::string &filename, bool simplex=false,
+ bool writeNodalValues=true) const
{
-
-// const bool simplex=false;
-
- FILE *f = fopen (filename.c_str(), "w");
- fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
+ FILE *f = fopen(filename.c_str(), "w");
+ fprintf(f, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
{
- fprintf(f,"$Nodes\n%d\n",_nodalValues.size());
+ fprintf(f, "$Nodes\n%d\n", (int)_nodalValues.size());
typename std::map<int, scalar>::const_iterator it = _nodalValues.begin();
- for ( ; it!=_nodalValues.end();++it){
+ for ( ; it != _nodalValues.end(); ++it){
SPoint3 p = coordinates_of_node(it->first);
- fprintf(f,"%d %g %g %g\n",it->first,p.x(),p.y(),p.z());
+ fprintf(f, "%d %g %g %g\n", it->first, p.x(), p.y(), p.z());
}
- fprintf(f,"$EndNodes\n");
+ fprintf(f, "$EndNodes\n");
}
{
- int coeff=1;
- if(simplex) coeff=6;
- fprintf(f,"$Elements\n%d\n",coeff*_active.size());
+ 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){
+ for ( ; it != _active.end(); ++it){
+ int i, j, k;
+ element_ijk(*it, i, j, k);
if(!simplex){
- fprintf(f,"%d 5 3 1 1 1",*it);
- int i,j,k;
- element_ijk(*it,i,j,k);
- 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,"\n");
- }else{
- int i,j,k;
- element_ijk(*it,i,j,k);
-
- //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,"\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,"\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,"\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,"\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,"\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,"\n");
+ 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, "\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, "\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, "\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, "\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, "\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, "\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, "\n");
}
}
fprintf(f,"$EndElements\n");
}
- if(writeLS){
- fprintf(f,"$NodeData\n1\n\"distance\"\n 1\n 0.0\n3\n0\n 1\n %d\n",_nodalValues.size());
+ if(writeNodalValues){
+ fprintf(f,"$NodeData\n1\n\"distance\"\n 1\n 0.0\n3\n0\n 1\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);
- // fprintf(f,"%d %g\n",it->first,distance(p.x(),p.y(),p.z()));
- fprintf(f,"%d %g\n",it->first,it->second);
- }
- fprintf(f,"$EndNodeData\n");
+ for( ; it != _nodalValues.end(); ++it)
+ fprintf(f, "%d %g\n", it->first, it->second);
+ fprintf(f, "$EndNodeData\n");
}
- fclose (f);
+ fclose(f);
}
-
- void writeLSOnly (const std::string &filename) const
+ void writeNodalValues(const std::string &filename) const
{
- FILE *f = fopen (filename.c_str(), "w");
-
- fprintf(f,"%d\n",_nodalValues.size());
+ FILE *f = fopen(filename.c_str(), "w");
+ fprintf(f, "%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);
- fprintf(f,"%d %g\n",it->first,it->second);
- }
-
- fclose (f);
-
+ for( ; it != _nodalValues.end(); ++it)
+ fprintf(f, "%d %g\n", it->first, it->second);
+ fclose(f);
}
-
};
-
-
#endif
diff --git a/utils/api_demos/mainCartesian.cpp b/utils/api_demos/mainCartesian.cpp
index 3d4b391847..de3f52327d 100644
--- a/utils/api_demos/mainCartesian.cpp
+++ b/utils/api_demos/mainCartesian.cpp
@@ -1,4 +1,3 @@
-#if 1
#include "Gmsh.h"
#include "GModel.h"
#include "MVertex.h"
@@ -6,228 +5,145 @@
#include "MTriangle.h"
#include "SOrientedBoundingBox.h"
#include "Numeric.h"
-#else
-#include <gmsh/Gmsh.h>
-#include <gmsh/GModel.h>
-#include <gmsh/MVertex.h>
-#include <gmsh/cartesian.h>
-#include <gmsh/MTriangle.h>
-#include <gmsh/SOrientedBoundingBox.h>
-#include <gmsh/Numeric.h>
-#endif
-void insertBoxes ( double x, double y, double z, double EP, cartesianBox<double> &box){
-
- int id1 = box.index_of_element(x-EP,y-EP,z-EP);
- int id2 = box.index_of_element(x+EP,y+EP,z+EP);
- int i1,j1,k1;
- box.element_ijk(id1,i1,j1,k1);
- int i2,j2,k2;
- box.element_ijk(id2,i2,j2,k2);
- // if (y > 0.9)printf("%g %g %g %3d %3d %3d %3d %3d %3d\n",x,y,z,i1,j1,k1,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);
- }
- }
- }
-}
-
-void test()
+void insertBoxes(double x, double y, double z, double thickness, cartesianBox<double> &box)
{
- printf("enter z coordinate : ");
- double x,y,z=-1;
- scanf ("%lf",&z);
- printf(" z = %g\n",z);
- std::vector<SPoint3> v, pts;
- std::vector<SPoint3> cv;
- std::vector<double> d;
-
-
- for (int i=0;i<100;i++){
- for (int j=0;j<100;j++){
- x = -.5+(double) i/50.;
- y = -.5+(double) j/50.;
- v.push_back(SPoint3(4.3,x,y));
+ 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 i1, j1 ,k1;
+ box.element_ijk(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);
+ }
}
}
- signedDistancesPointsTriangle (d,cv,v,SPoint3(4,0,0),SPoint3(4,1,0),SPoint3(4,0,1));
-
- FILE *f = fopen ("chitte.pos","w");
- fprintf(f,"View \"\"{\n");
- for( int i=0;i<v.size();i++){
- fprintf(f,"SP(%g,%g,%g){%g};\n",v[i].x(),v[i].y(),v[i].z(),d[i]);
- }
- fprintf(f,"};\n");
- fclose(f);
}
int main (int argc,char *argv[])
{
- // test();
- // return 1;
-
- // ne pas oublier les conditions limites
-
- GmshInitialize();
-
- if (argc != 6){
- printf("usage : mainCartesian meshFile thickness NPointsX SAMPLING(smaller than thickness) filterCells\n");
+ if(argc < 4){
+ printf("Usage: %s file thickness npointsx [sampling] [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");
return -1;
}
- const int FACT=atoi(argv[3]);
- const int FILTER=atoi(argv[5]);
-
- double EP = atof(argv[2]);
-
- double SAMPLING;
- if (argc == 5)SAMPLING = atof(argv[4]);
- else SAMPLING = EP/3.0;
-
+ GmshInitialize();
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;
+
GModel *gm = GModel::current();
- printf("mesh read\n");
-
- std::vector<SPoint3> POINTS, dummy;
- std::vector<SVector3> NORMALS;
+ std::vector<SPoint3> points;
+ std::vector<SVector3> normals;
- for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
- (*fit)->fillPointCloud(SAMPLING,&POINTS,&NORMALS);
- }
+ 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());
- // printf("OBB with %d points\n",POINTS.size());
- // SOrientedBoundingBox *obb = SOrientedBoundingBox::buildOBB (POINTS);
- // printf("OBB done\n");
+ // printf("OBB with %d points\n", points.size());
+ // SOrientedBoundingBox *obb = SOrientedBoundingBox::buildOBB(points);
+ // printf("OBB done\n");
SBoundingBox3d bb;
- for(int i=0;i<POINTS.size();i++){
- bb += POINTS[i];
- }
-
- bb.scale(1.2,1.2,1.2);
-
-
+ for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
+ bb.scale(1.2, 1.2, 1.2);
SVector3 range = bb.max() - bb.min();
- int NX =FACT;
- int NY =FACT*range.y()/range.x();
- int NZ =FACT*range.z()/range.x();
-
- 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);
-
-
- 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);
-
- box.points() = POINTS;
- box.normals() = NORMALS;
-
- printf("%d points in the cloud \n", box.points().size());
+ int NX = npointsx;
+ int NY = NX * range.y() / range.x();
+ int NZ = NX * range.z() / range.x();
-
- // FILE *f = fopen ("chitte.pos","w");
- // fprintf(f,"View \"\"{\n");
+ 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);
- for (int i=0; i< POINTS.size();i++){
- insertBoxes ( POINTS [i].x(), POINTS [i].y(), POINTS [i].z(), EP, box);
- // fprintf(f,"SP(%g,%g,%g){%g};\n",POINTS[i].x(),POINTS[i].y(),POINTS[i].z(),0.0);
- }
- // fprintf(f,"};\n");
- // fclose(f);
+ 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);
+ for (unsigned int i = 0; i < points.size();i++)
+ insertBoxes(points[i].x(), points[i].y(), points[i].z(), thickness, box);
printf("insertion done\n");
box.create_nodes();
-
- std::map<int,double>::const_iterator it =
- box.begin();
+ 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));
+ 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",NODES.size());
- int III = 0;
+ printf("%d nodes in the cartesian mesh\n", (int)NODES.size());
+
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
- // printf("model face %d : %d triangles\n",(*fit)->tag(),(*fit)->stl_triangles.size()/3);
- for (int i=0;i<(*fit)->stl_triangles.size(); i+=3){
+ 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];
+ 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])*.3333333333;
+ 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));
-
- // printf("N1 %g %g %g N2 %g %g %g -- %g %g %g -- %g %g %g -- %g %g %g\n",N.x(),N.y(),N.z(),NN.x(),NN.y(),NN.z(),
- // p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),p3.x(),p3.y(),p3.z());
+ 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);
+ signedDistancesPointsTriangle(localdist, CNODES, NODES, P1, P2, P3);
else
- signedDistancesPointsTriangle (localdist,CNODES,NODES,P2,P1,P3);
-
- if(1){
- if (dist.empty()) dist=localdist;
- else
- for (int j=0;j<localdist.size();j++)
- dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
- }
+ 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];
}
- III++;
}
-
- for (int j=0;j<dist.size();j++){
- // printf("d(%d)=%g\n",indices[j],dist[j]);
- box.setValue(indices[j],dist[j]);
- }
-
+ for (unsigned int j = 0; j < dist.size(); j++)
+ box.setValue(indices[j], dist[j]);
-#if 1
- if(FILTER){
- int nbErased=0;
+ if(filter){
+ int nbErased = 0;
//Coup de menage avant d'exporter le maillage
- for( cartesianBox<double>::boxIter it=box.activeBoxBegin(); it!=box.activeBoxEnd();){
-
+ for(cartesianBox<double>::boxIter it = box.activeBoxBegin(); it != box.activeBoxEnd();){
std::vector<double> ls_vals;
box.getNodalValues(*it, ls_vals);
-
- double lsmax= *std::max_element(ls_vals.begin(), ls_vals.end());
+ double lsmax = *std::max_element(ls_vals.begin(), ls_vals.end());
double lsmin = *std::min_element(ls_vals.begin(), ls_vals.end());
- double change_sign = lsmax*lsmin ;
-
- // std::cout<<"active cell is "<<*it<<" and change sign is "<<change_sign<<" lsmax is "<<lsmax<<std::endl;
+ double change_sign = lsmax * lsmin ;
double epsilon = 1.e-10;
- if(change_sign>0 && lsmax < -epsilon) {box.erase(*(it++)); nbErased++;}
+ if(change_sign > 0 && lsmax < -epsilon) {
+ box.erase(*(it++));
+ nbErased++;
+ }
else ++it;
}
-
- std::cout<<"Number of erased cells after filtering : "<<nbErased<<std::endl;
+ std::cout << "Number of erased cells after filtering : " << nbErased << std::endl;
}
-#endif
printf("nodes created\n");
- box.writeMSH("yeah.msh",true);
+ box.writeMSH("yeah.msh", true);
- box.writeMSH("youhou.msh",true,false);
- box.writeLSOnly("youhou.pos");
+ box.writeMSH("youhou.msh", true, false);
+ box.writeNodalValues("youhou.pos");
printf("mesh written\n");
GmshFinalize();
}
--
GitLab