Skip to content
Snippets Groups Projects
Commit eed986f5 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

No commit message 

No commit message
parent fab5a46d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/cartesian.h
+ 58
− 60
View file @ eed986f5
......@@ -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());
......
doc/CREDITS.txt
+ 1
− 1
View file @ eed986f5
Gmsh is copyright (C) 1997-2009
Gmsh is copyright (C) 1997-2010
Christophe Geuzaine
<cgeuzaine at ulg.ac.be>
......
utils/api_demos/mainCartesian.cpp
+ 43
− 40
View file @ eed986f5
......@@ -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");
......
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