From eed986f556c08cfe5b2267c5783a5151fb4471c7 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 21 Jul 2010 09:59:54 +0000
Subject: [PATCH]
---
Numeric/cartesian.h | 118 +++++++++++++++---------------
doc/CREDITS.txt | 2 +-
utils/api_demos/mainCartesian.cpp | 83 +++++++++++----------
3 files changed, 102 insertions(+), 101 deletions(-)
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index bd6f0d238e..68b82937b7 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -8,50 +8,57 @@
#include "SVector3.h"
#include "SPoint3.h"
-// This is a cartesian mesh that encompasses an oriented box with NXI
-// * NETA * NZETA hexaderal cells, with values stored at vertices
-
-/*
- 2-D example:
-
- j
- +---+---+---+---+---+---+
- | | | |ij | | |
- i +---+---+---+---+---+---+
- | | | | | | |
- +---+---+---+---+---+---+
-
- 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
-*/
-
+// A cartesian grid that encompasses an oriented 3-D box, with values
+// stored at vertices:
+//
+// j
+// +---+---+---+---+---+---+
+// | | | | | | |
+// i +---+---+-(i,j)-+---+---+
+// | | | | | | |
+// +---+---+---+---+---+---+
+//
+// Nodal values and active hexahedral cells are stored and
+// referenced by a linear index (i + N * j)
+//
+// The (i,j) cell has nodes (i,j), (i+1,j), (i+1,j+1) and (i,j+1)
template <class scalar>
class cartesianBox {
private:
+ // number of subdivisions along the xi-, eta- and zeta-axis
int _Nxi, _Neta, _Nzeta;
+ // origin of the grid and spacing along xi, eta and zeta
double _X, _Y, _Z, _dxi, _deta, _dzeta;
+ // xi-, eta- and zeta-axis directions
SVector3 _xiAxis, _etaAxis, _zetaAxis;
+ // set of active cells; the value stored for cell (i,j,k) is the
+ // linear index (i + _Nxi * j + _Nxi *_Neta * k)
std::set<int> _activeCells;
+ // map of stored nodal values, index by the linear index
+ // (i + (_Nxi+1) * j + (_Nxi+1) * (_Neta+1) * k)
typename std::map<int, scalar> _nodalValues;
+ // mapping from linear node index to unique node tags
+ std::map<int, int> _nodeTags;
+ // level of the box (coarset box has highest level; finest box has
+ // level==1)
int _level;
+ // pointer to a finer (refined by 2) level box (if any)
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 level=1)
+ const SVector3 &dxi, const SVector3 &deta, const SVector3 &dzeta,
+ 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), _level(level), _childBox(0)
+ _dxi(norm(dxi)), _deta(norm(deta)), _dzeta(norm(dzeta)),
+ _xiAxis(dxi), _etaAxis(deta), _zetaAxis(dzeta),
+ _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,
+ _childBox = new cartesianBox<scalar>(X, Y, Z, dxi, deta, dzeta,
+ 2 * Nxi, 2 * Neta, 2 * Nzeta,
level - 1);
}
int getNxi(){ return _Nxi; }
@@ -66,38 +73,25 @@ class cartesianBox {
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)
+ void getCellValues(int t, std::vector<scalar> &values)
{
- // perhaps not optimal (lot of searches...)
int i, j, k;
getCellIJK(t, i, j, k);
- typename std::map<int, scalar>::iterator itNode;
-
- itNode = _nodalValues.find(getNodeIndex(i, j, k));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i+1,j,k));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i+1,j+1,k));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i,j+1,k));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i,j,k+1));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i+1,j,k+1));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i+1,j+1,k+1));
- ls_values.push_back(itNode->second);
-
- itNode = _nodalValues.find(getNodeIndex(i,j+1,k+1));
- ls_values.push_back(itNode->second);
+ for(int I = 0; I < 2; I++)
+ for(int J = 0; J < 2; J++)
+ for(int K = 0; K < 2; K++){
+ typename std::map<int, scalar>::iterator it =
+ _nodalValues.find(getNodeIndex(i + I, j + J, k + K));
+ if(it != _nodalValues.end())
+ values.push_back(it->second);
+ else{
+ Msg::Error("Could not find value i,j,k=%d,%d,%d for cell %d\n",
+ i + I, j + J, k + K, t);
+ values.push_back(0.);
+ }
+ }
}
- inline int getCellContainingPoint(double x, double y, double z) const
+ int getCellContainingPoint(double x, double y, double z) const
{
// P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
// DP = P-P_0 * _vdx = xi
@@ -113,7 +107,7 @@ class cartesianBox {
if (k < 0) k = 0; if (k >= _Nzeta) k = _Nzeta - 1;
return getCellIndex(i, j, k);
}
- inline SPoint3 getNodeCoordinates(const int &t) const
+ SPoint3 getNodeCoordinates(const int &t) const
{
int i, j, k;
getNodeIJK(t, i, j, k);
@@ -130,27 +124,27 @@ class cartesianBox {
{
return (_activeCells.find(t) != _activeCells.end());
}
- inline int getCellIndex(int i, int j, int k) const
+ int getCellIndex(int i, int j, int k) const
{
return i + _Nxi * j + _Nxi *_Neta * k;
}
- inline int getNodeIndex(int i, int j, int k) const
+ int getNodeIndex(int i, int j, int k) const
{
return i + (_Nxi+1) * j + (_Nxi+1) * (_Neta+1) * k;
}
- inline void getCellIJK(int index, int &i, int &j, int &k) const
+ 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 getNodeIJK(int index, int &i, int &j, int &k) const
+ 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 createNodalValues()
+ void createNodalValues()
{
std::set<int>::const_iterator it = _activeCells.begin();
for( ; it != _activeCells.end() ; ++it){
@@ -168,6 +162,10 @@ class cartesianBox {
bool writeNodalValues=true) const
{
FILE *f = fopen(filename.c_str(), "w");
+ if(!f){
+ Msg::Error("Could not open file '%s'", filename.c_str());
+ return;
+ }
fprintf(f, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
{
fprintf(f, "$Nodes\n%d\n", (int)_nodalValues.size());
diff --git a/doc/CREDITS.txt b/doc/CREDITS.txt
index 3d30a59d72..41e60ebb27 100644
--- a/doc/CREDITS.txt
+++ b/doc/CREDITS.txt
@@ -1,4 +1,4 @@
- Gmsh is copyright (C) 1997-2009
+ Gmsh is copyright (C) 1997-2010
Christophe Geuzaine
<cgeuzaine at ulg.ac.be>
diff --git a/utils/api_demos/mainCartesian.cpp b/utils/api_demos/mainCartesian.cpp
index 052ce22e35..386fa294fe 100644
--- a/utils/api_demos/mainCartesian.cpp
+++ b/utils/api_demos/mainCartesian.cpp
@@ -77,8 +77,7 @@ void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
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));
+ points.push_back(SPoint3(p.x(), p.y(), p.z()));
}
}
@@ -143,10 +142,32 @@ void removeParentCellsWithChildren(cartesianBox<double> *box)
removeParentCellsWithChildren(box->getChildBox());
}
+void removeOutsideCells(cartesianBox<double> *box)
+{
+ if(!box) return;
+ int nbErased = 0;
+ for(cartesianBox<double>::cellIter it = box->activeCellsBegin();
+ it != box->activeCellsEnd();){
+ std::vector<double> vals;
+ box->getCellValues(*it, vals);
+ double lsmax = *std::max_element(vals.begin(), vals.end());
+ double lsmin = *std::min_element(vals.begin(), vals.end());
+ double change_sign = lsmax * lsmin;
+ double epsilon = 1.e-10;
+ if(change_sign > 0 && lsmax < -epsilon) {
+ box->eraseActiveCell(*(it++));
+ nbErased++;
+ }
+ else ++it;
+ }
+ Msg::Info(" number of cells erased after filtering: %d", nbErased);
+ removeOutsideCells(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("Usage: %s file lx ly lz rmax [levels=1]\n", argv[0]);
printf("where\n");
printf(" 'file' contains a CAD model\n");
printf(" 'lx', 'ly' and 'lz' are the sizes of the elements along the"
@@ -154,8 +175,6 @@ int main(int argc,char *argv[])
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;
}
@@ -165,10 +184,9 @@ int main(int argc,char *argv[])
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.;
+ double sampling = std::min(rmax, 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.;
@@ -176,17 +194,17 @@ int main(int argc,char *argv[])
GModel *gm = GModel::current();
std::vector<SPoint3> points;
- Msg::Info("Filling point clould on surfaces");
+ Msg::Info("Filling point cloud on surfaces");
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
(*fit)->fillPointCloud(sampling, &points);
Msg::Info(" %d points in the surface cloud", (int)points.size());
std::vector<SPoint3> edgePoints;
if(levels > 1){
- Msg::Info("Filling point clould on curves");
+ Msg::Info("Filling point cloud 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());
+ Msg::Info(" %d points in the curve cloud", (int)edgePoints.size());
}
SBoundingBox3d bb;
@@ -201,10 +219,11 @@ int main(int argc,char *argv[])
if(NY < 2) NY = 2;
if(NZ < 2) NZ = 2;
- Msg::Info(" bb Min= %g %g %g -- bb Max= %g %g %g -- NX %d NY %d NZ %d",
+ Msg::Info(" bounding box min: %g %g %g -- max: %g %g %g",
bb.min().x(), bb.min().y(), bb.min().z(),
- bb.max().x(), bb.max().y(), bb.max().z(), NX, NY, NZ);
-
+ bb.max().x(), bb.max().y(), bb.max().z());
+ Msg::Info(" Nx=%d Ny=%d Nz=%d", 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),
@@ -222,43 +241,27 @@ int main(int argc,char *argv[])
parent = child;
}
- 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
+ // which there is no parent neighbor; then remove parent cells that
+ // have children
+ Msg::Info("Removing cells to match X-FEM mesh topology constraints");
removeOrphanChildCells(&box);
- // remove parent cells that have children
removeParentCellsWithChildren(&box);
- // TODO:
- // * remove child nodes which exist in parent grid
- // * offset of node numbers depending on box level to write a single mesh
-
+ // we generate duplicate nodes at this point so we can easily access
+ // cell values at each level; we will clean up by renumbering after
+ // filtering
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>::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());
- double lsmin = *std::min_element(ls_vals.begin(), ls_vals.end());
- double change_sign = lsmax * lsmin ;
- double epsilon = 1.e-10;
- if(change_sign > 0 && lsmax < -epsilon) {
- box.eraseActiveCell(*(it++));
- nbErased++;
- }
- else ++it;
- }
- Msg::Info(" number of erased cells after filtering: %d", nbErased);
- }
+ Msg::Info("Removing cells outside the structure");
+ removeOutsideCells(&box);
+
+ Msg::Info("Renumbering mesh vertices across levels");
+
Msg::Info("Writing results to disk");
box.writeMSH("yeah.msh");
--
GitLab