diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index b9967b6b092c36d99c32b0295558328dbbf46d5b..f849a8b18788c4116ae815c57400397b2a3d8857 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -85,11 +85,12 @@ void PrintUsage(const char *name)
Msg::Direct(" -check Perform various consistency checks on mesh");
#if defined(HAVE_FLTK)
Msg::Direct("Post-processing options:");
- Msg::Direct(" -noview Hide all views on startup");
Msg::Direct(" -link int Select link mode between views (0, 1, 2, 3, 4)");
Msg::Direct(" -combine Combine views having identical names into multi-time-step views");
Msg::Direct("Display options:");
Msg::Direct(" -nodb Disable double buffering");
+ Msg::Direct(" -noview Hide all views on startup");
+ Msg::Direct(" -nomesh Hide all meshes on startup");
Msg::Direct(" -fontsize int Specify the font size for the GUI");
Msg::Direct(" -theme string Specify FLTK GUI theme");
Msg::Direct(" -display string Specify display");
@@ -580,6 +581,15 @@ void GetOptions(int argc, char *argv[])
opt_view_visible(0, GMSH_SET, 0);
i++;
}
+ else if(!strcmp(argv[i] + 1, "nomesh")) {
+ opt_mesh_points(0, GMSH_SET, 0.);
+ opt_mesh_lines(0, GMSH_SET, 0.);
+ opt_mesh_surfaces_edges(0, GMSH_SET, 0.);
+ opt_mesh_surfaces_faces(0, GMSH_SET, 0.);
+ opt_mesh_volumes_edges(0, GMSH_SET, 0.);
+ opt_mesh_volumes_faces(0, GMSH_SET, 0.);
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "link")) {
i++;
if(argv[i])
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index c53b3f4ce755602087c7323b887ed1df726489a4..ec953370ddfbe53a44266035258c38de870e5d4d 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1182,16 +1182,6 @@ void GFace::moveToValidRange(SPoint2 &pt) const
}
}
-void GFace::addTriangle(MTriangle *t) {
- triangles.push_back(t);
-}
-void GFace::addQuadrangle(MQuadrangle *q) {
- quadrangles.push_back(q);
-}
-void GFace::addPolygon(MPolygon *p) {
- polygons.push_back(p);
-}
-
#include "Bindings.h"
void GFace::registerBindings(binding *b)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 9f71f5981bc8b25c8af1df17f67a699954bb70ff..00aee5a100a99577dff8b6b3ceae1d0dc568a4bf 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -299,9 +299,9 @@ class GFace : public GEntity
std::vector<MQuadrangle*> quadrangles;
std::vector<MPolygon*> polygons;
- void addTriangle(MTriangle *t);
- void addQuadrangle(MQuadrangle *q);
- void addPolygon(MPolygon *p);
+ void addTriangle(MTriangle *t){ triangles.push_back(t); }
+ void addQuadrangle(MQuadrangle *q){ quadrangles.push_back(q); }
+ void addPolygon(MPolygon *p){ polygons.push_back(p); }
// an array with additional vertices that are supposed to exist
// in the final mesh of the model face. This can be used for
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 843b2a0e22ef424e2265c0e49bb889121d623f32..e59e05e5e64192f4cdf849cc39753f8898468bdd 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -5,6 +5,7 @@
#ifndef _DIVIDE_AND_CONQUER_H_
#define _DIVIDE_AND_CONQUER_H_
+
#include <vector>
#include <algorithm>
#include "fullMatrix.h"
diff --git a/Numeric/Gauss.h b/Numeric/Gauss.h
index 6be921404359def80fcfc4479cfc6ef8ebcefd2e..06ee437403d7b3a14af093f00022a090bc32ca94 100644
--- a/Numeric/Gauss.h
+++ b/Numeric/Gauss.h
@@ -5,8 +5,8 @@
#ifndef _GAUSS_H_
#define _GAUSS_H_
-#include "fullMatrix.h"
+#include "fullMatrix.h"
struct IntPt{
double pt[3];
@@ -36,7 +36,9 @@ IntPt *getGQPriPts(int order);
int getNGQHPts(int order);
IntPt *getGQHPts(int order);
-//For now this class is only for bindings but maybe the interface is cleaner (it does not add new types) and it can replace the other interface
+//For now this class is only for bindings but maybe the interface is
+//cleaner (it does not add new types) and it can replace the other
+//interface
class gaussIntegration {
public:
static void getTriangle(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index 68b82937b739fc9535a0a92cc8b658f2d297b881..09c6bc6904a051588221bde2599d9077615a1e8f 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _CARTESIAN_H_
#define _CARTESIAN_H_
@@ -7,6 +12,7 @@
#include <stdio.h>
#include "SVector3.h"
#include "SPoint3.h"
+#include "GmshMessage.h"
// A cartesian grid that encompasses an oriented 3-D box, with values
// stored at vertices:
@@ -24,7 +30,7 @@
// 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:
+ 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
@@ -34,16 +40,91 @@ class cartesianBox {
// 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;
+ // map of stored nodal values, index by the linear index (i +
+ // (_Nxi+1) * j + (_Nxi+1) * (_Neta+1) * k). Along with the value is
+ // stored a node tag (used for global numbering of the nodes across
+ // the grid levels)
+ typename std::map<int, std::pair<scalar, int> > _nodalValues;
// 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;
+ int _getNumNodes()
+ {
+ int n = 0;
+ for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); it++)
+ if(it->second.second > 0) n++;
+ if(_childBox) n += _childBox->_getNumNodes();
+ return n;
+ }
+ void _printNodes(FILE *f)
+ {
+ for (valIter it = _nodalValues.begin(); it != _nodalValues.end(); ++it){
+ if(it->second.second > 0){
+ SPoint3 p = getNodeCoordinates(it->first);
+ fprintf(f, "%d %g %g %g\n", it->second.second, p.x(), p.y(), p.z());
+ }
+ }
+ if(_childBox) _childBox->_printNodes(f);
+ }
+ int _getNumElements(bool simplex)
+ {
+ int coeff = simplex ? 6 : 1;
+ int n = _activeCells.size() * coeff;
+ if(_childBox) n += _childBox->_getNumElements(simplex);
+ return n;
+ }
+ void _printElements(FILE *f, bool simplex, int startingNum=1)
+ {
+ int num = startingNum;
+ for(cellIter it = _activeCells.begin(); it != _activeCells.end(); ++it){
+ int i, j, k;
+ getCellIJK(*it, i, j, k);
+ if(!simplex){
+ fprintf(f, "%d 5 3 1 1 1", num++);
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i, j, k))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i + 1, j, k))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i + 1, j + 1, k))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i, j + 1, k))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i, j, k + 1))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i + 1, j, k + 1))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i + 1, j + 1, k + 1))));
+ fprintf(f, " %d", std::abs(getNodeTag(getNodeIndex(i, j + 1, k + 1))));
+ fprintf(f, "\n");
+ }
+ else{
+ int idx[6][4] = {
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i, j + 1, k + 1),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j + 1, k + 1)},
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j + 1, k + 1),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j + 1, k)},
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i, j, k + 1),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j + 1, k + 1)},
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j + 1, k),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i + 1, j, k)},
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i + 1, j, k),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j, k)},
+ {getNodeIndex(i, j + 1, k), getNodeIndex(i, j, k),
+ getNodeIndex(i + 1, j, k + 1), getNodeIndex(i, j, k + 1)}
+ };
+ for(int ii = 0; ii < 6; ii++){
+ fprintf(f, "%d 4 3 1 1 1 %d %d %d %d\n", num++,
+ std::abs(getNodeTag(idx[ii][0])), std::abs(getNodeTag(idx[ii][1])),
+ std::abs(getNodeTag(idx[ii][2])), std::abs(getNodeTag(idx[ii][3])));
+ }
+ }
+ }
+ if(_childBox) _childBox->_printElements(f, simplex, num);
+ }
+ void _printValues(FILE *f)
+ {
+ for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); ++it){
+ if(it->second.second > 0)
+ fprintf(f, "%d %g\n", it->second.second, it->second.first);
+ }
+ if(_childBox) _childBox->_printValues(f);
+ }
public:
cartesianBox(double X, double Y, double Z,
const SVector3 &dxi, const SVector3 &deta, const SVector3 &dzeta,
@@ -69,21 +150,20 @@ class cartesianBox {
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;
+ typedef typename std::map<int, std::pair<scalar, int> >::iterator valIter;
valIter nodalValuesBegin(){ return _nodalValues.begin(); }
valIter nodalValuesEnd(){ return _nodalValues.end(); }
- void setNodalValue(int i, scalar s){ _nodalValues[i] = s; }
- void getCellValues(int t, std::vector<scalar> &values)
+ void setNodalValue(int i, scalar s){ _nodalValues[i].first = s; }
+ void getNodalValuesForCell(int t, std::vector<scalar> &values)
{
int 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++){
- typename std::map<int, scalar>::iterator it =
- _nodalValues.find(getNodeIndex(i + I, j + J, k + K));
+ valIter it = _nodalValues.find(getNodeIndex(i + I, j + J, k + K));
if(it != _nodalValues.end())
- values.push_back(it->second);
+ values.push_back(it->second.first);
else{
Msg::Error("Could not find value i,j,k=%d,%d,%d for cell %d\n",
i + I, j + J, k + K, t);
@@ -107,7 +187,7 @@ class cartesianBox {
if (k < 0) k = 0; if (k >= _Nzeta) k = _Nzeta - 1;
return getCellIndex(i, j, k);
}
- SPoint3 getNodeCoordinates(const int &t) const
+ SPoint3 getNodeCoordinates(int t) const
{
int i, j, k;
getNodeIJK(t, i, j, k);
@@ -117,10 +197,9 @@ class cartesianBox {
SVector3 D = xi * _xiAxis + eta * _etaAxis + zeta * _zetaAxis;
return SPoint3(_X + D.x(), _Y + D.y(), _Z + D.z());
}
- 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)
+ void insertActiveCell(int t){ _activeCells.insert(t); }
+ void eraseActiveCell(int t){ _activeCells.erase(t); }
+ bool activeCellExists(int t)
{
return (_activeCells.find(t) != _activeCells.end());
}
@@ -132,6 +211,12 @@ class cartesianBox {
{
return i + (_Nxi+1) * j + (_Nxi+1) * (_Neta+1) * k;
}
+ int getNodeTag(int index)
+ {
+ valIter it = _nodalValues.find(index);
+ if(it != _nodalValues.end()) return it->second.second;
+ else return 0;
+ }
void getCellIJK(int index, int &i, int &j, int &k) const
{
k = index / (_Nxi * _Neta);
@@ -146,123 +231,61 @@ class cartesianBox {
}
void createNodalValues()
{
- std::set<int>::const_iterator it = _activeCells.begin();
- for( ; it != _activeCells.end() ; ++it){
+ for(cellIter it = _activeCells.begin(); it != _activeCells.end(); ++it){
const int &t = *it;
int 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[getNodeIndex(i + I, j + J, k + K)] = 0.;
+ _nodalValues[getNodeIndex(i + I, j + J, k + K)] =
+ std::pair<scalar, int>(0., 0);
}
if(_childBox) _childBox->createNodalValues();
}
- void writeMSH(const std::string &filename, bool simplex=false,
- bool writeNodalValues=true) const
+ void renumberNodes(int startingNum=1, cartesianBox<scalar> *parent=0)
+ {
+ int num = startingNum;
+ for(valIter it = _nodalValues.begin(); it != _nodalValues.end(); it++){
+ int i, j, k;
+ getNodeIJK(it->first, i, j, k);
+ if(!parent || i % 2 || j % 2 || k % 2)
+ it->second.second = num++;
+ else{
+ int tag = parent->getNodeTag(parent->getNodeIndex(i / 2, j / 2, k / 2));
+ if(!tag) // FIXME! not sure why this can happen, but it does (bug?)
+ it->second.second = num++;
+ else // the node exists in the coarset grid: store it with negative sign
+ it->second.second = -std::abs(tag);
+ }
+ }
+ if(_childBox) _childBox->renumberNodes(num, this);
+ }
+ void writeMSH(const std::string &fileName, bool simplex=false,
+ bool writeNodalValues=true)
{
- FILE *f = fopen(filename.c_str(), "w");
+ FILE *f = fopen(fileName.c_str(), "w");
if(!f){
- Msg::Error("Could not open file '%s'", filename.c_str());
+ Msg::Error("Could not open file '%s'", fileName.c_str());
return;
}
+ int numNodes = _getNumNodes(), numElements = _getNumElements(simplex);
+ Msg::Info("Writing '%s' (%d nodes, %d elements)", fileName.c_str(),
+ numNodes, numElements);
fprintf(f, "$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
- {
- 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 = 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)_activeCells.size());
- std::set<int>::const_iterator it = _activeCells.begin();
- for ( ; it != _activeCells.end(); ++it){
- int i, j, k;
- getCellIJK(*it, i, j, k);
- if(!simplex){
- fprintf(f, "%d 5 3 1 1 1", *it);
- 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", 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", 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", 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", 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", 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", 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);
- }
+ fprintf(f, "$Nodes\n%d\n", numNodes);
+ _printNodes(f);
+ fprintf(f, "$EndNodes\n");
+ fprintf(f,"$Elements\n%d\n", numElements);
+ _printElements(f, simplex);
+ fprintf(f,"$EndElements\n");
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)
- fprintf(f, "%d %g\n", it->first, it->second);
+ fprintf(f,"$NodeData\n1\n\"distance\"\n1\n0.0\n3\n0\n1\n%d\n", numNodes);
+ _printValues(f);
fprintf(f, "$EndNodeData\n");
}
fclose(f);
}
- void writeNodalValues(const std::string &filename) const
- {
- 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)
- fprintf(f, "%d %g\n", it->first, it->second);
- fclose(f);
- }
};
#endif
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index d65901bb94309329bbb44018abcc0be61aebb117..00698341af61c76526e7b2b702f076ca0d54919e 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -888,8 +888,6 @@ Perform various consistency checks on mesh
@noindent Post-processing options:
@ftable @code
-@item -noview
-Hide all views on startup
@item -link int
Select link mode between views (0, 1, 2, 3, 4)
@item -combine
@@ -903,6 +901,10 @@ Combine views having identical names into multi-time-step views
@ftable @code
@item -nodb
Disable double buffering
+@item -noview
+Hide all views on startup
+@item -nomesh
+Hide all meshes on startup
@item -fontsize int
Specify the font size for the GUI
@item -theme string
diff --git a/utils/api_demos/mainCartesian.cpp b/utils/api_demos/mainCartesian.cpp
index 386fa294fe046924fbc43ac9047106b40faba5bc..cc8a8d0c8f8bce00260832c0f9f51ac2b240fec8 100644
--- a/utils/api_demos/mainCartesian.cpp
+++ b/utils/api_demos/mainCartesian.cpp
@@ -6,7 +6,6 @@
#include "SOrientedBoundingBox.h"
#include "Numeric.h"
#include "GmshMessage.h"
-#include "OS.h"
void insertActiveCells(double x, double y, double z, double rmax,
cartesianBox<double> &box)
@@ -32,7 +31,7 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
nodes.push_back(box.getNodeCoordinates(it->first));
indices.push_back(it->first);
}
- Msg::Info(" %d nodes in the grid", (int)nodes.size());
+ Msg::Info(" %d nodes in the grid at level %d", (int)nodes.size(), box.getLevel());
std::vector<double> dist, localdist;
std::vector<SPoint3> dummy;
for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
@@ -81,7 +80,7 @@ void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
}
}
-int removeOrphanChildCells(cartesianBox<double> *parent)
+int removeBadChildCells(cartesianBox<double> *parent)
{
cartesianBox<double> *child = parent->getChildBox();
if(!child) return 0;
@@ -117,7 +116,7 @@ int removeOrphanChildCells(cartesianBox<double> *parent)
(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);
+ return removeBadChildCells(child);
}
void removeParentCellsWithChildren(cartesianBox<double> *box)
@@ -149,7 +148,7 @@ void removeOutsideCells(cartesianBox<double> *box)
for(cartesianBox<double>::cellIter it = box->activeCellsBegin();
it != box->activeCellsEnd();){
std::vector<double> vals;
- box->getCellValues(*it, vals);
+ box->getNodalValuesForCell(*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;
@@ -186,7 +185,7 @@ int main(int argc,char *argv[])
int levels = (argc > 6) ? atof(argv[6]) : 1;
// minimum distance between points in the cloud at the coarsest level
- double sampling = std::min(rmax, std::min(lx, std::min(ly, lz))) / 2.;
+ double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
// radius of the "tube" created around edges at the coarsest level
double rtube = std::max(lx, std::max(ly, lz)) * 2.;
@@ -194,23 +193,28 @@ int main(int argc,char *argv[])
GModel *gm = GModel::current();
std::vector<SPoint3> points;
- Msg::Info("Filling point cloud on surfaces");
+ Msg::Info("Filling coarse 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;
+ std::vector<SPoint3> refinePoints;
if(levels > 1){
- Msg::Info("Filling point cloud on curves");
+ double s = sampling / pow(2., levels - 1);
+ Msg::Info("Filling refined point cloud on curves and curved surfaces");
for (GModel::eiter eit = gm->firstEdge(); eit != gm->lastEdge(); eit++)
- fillPointCloud(*eit, sampling / pow(2., levels - 1), edgePoints);
- Msg::Info(" %d points in the curve cloud", (int)edgePoints.size());
+ fillPointCloud(*eit, s, refinePoints);
+ // FIXME: refine this by computing e.g. "mean" curvature
+ for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
+ if((*fit)->geomType() != GEntity::Plane)
+ (*fit)->fillPointCloud(2 * s, &refinePoints);
+ Msg::Info(" %d points in the refined cloud", (int)refinePoints.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);
+ for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
+ bb.scale(1.21, 1.21, 1.21);
SVector3 range = bb.max() - bb.min();
int NX = range.x() / lx;
int NY = range.y() / ly;
@@ -231,13 +235,16 @@ int main(int argc,char *argv[])
NX, NY, NZ, levels);
Msg::Info("Inserting active cells in the cartesian grid");
- for (unsigned int i = 0; i < points.size();i++)
+ Msg::Info(" level %d", box.getLevel());
+ for (unsigned int i = 0; i < points.size(); i++)
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);
+ Msg::Info(" level %d", child->getLevel());
+ for(unsigned int i = 0; i < refinePoints.size(); i++)
+ insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
+ rtube / pow(2., (levels - child->getLevel())), *child);
parent = child;
}
@@ -245,7 +252,7 @@ int main(int argc,char *argv[])
// 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);
+ removeBadChildCells(&box);
removeParentCellsWithChildren(&box);
// we generate duplicate nodes at this point so we can easily access
@@ -261,12 +268,10 @@ int main(int argc,char *argv[])
removeOutsideCells(&box);
Msg::Info("Renumbering mesh vertices across levels");
+ box.renumberNodes();
-
- Msg::Info("Writing results to disk");
- box.writeMSH("yeah.msh");
- box.writeMSH("youhou.msh", true, false);
- box.writeNodalValues("youhou.pos");
+ box.writeMSH("yeah.msh", true);
+
Msg::Info("Done!");
GmshFinalize();
}