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

cleanup

parent c686b931
No related branches found
No related tags found
No related merge requests found
...@@ -17,5 +17,5 @@ endif(GMSH_LIB MATCHES ".a") ...@@ -17,5 +17,5 @@ endif(GMSH_LIB MATCHES ".a")
include_directories(${GMSH_INC}) include_directories(${GMSH_INC})
add_executable(mainSimple mainSimple.cpp) add_executable(t0 t0.cpp)
target_link_libraries(mainSimple ${GMSH_LIB} ${LAPACK_LIB} ${BLAS_LIB}) target_link_libraries(t0 ${GMSH_LIB} ${LAPACK_LIB} ${BLAS_LIB})
//
// A simple example on how to build a GUI frontend to Gmsh using GLUT
// and libAntTweakBar
//
#include <cstring>
#if defined(__APPLE__)
# include <GLUT/glut.h>
#else
# include <GL/glut.h>
#endif
#include <AntTweakBar.h>
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "drawContext.h"
static drawContext *ctx = 0;
static mousePosition clickPos, prevPos;
static int specialKey = 0;
class drawContextTw : public drawContextGlobal{
public:
void draw(){ ctx->draw3d(); ctx->draw2d(); }
const char *getFontName(int index){ return "Helvetica"; }
int getFontSize(){ return 18; }
double getStringWidth(const char *str)
{
return glutBitmapLength(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)str);
}
int getStringHeight(){ return 18; }
int getStringDescent(){ return 6; }
void drawString(const char *str)
{
for (int i = 0; i < strlen(str); i++)
glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, str[i]);
}
};
// GLUT callbacks
void display()
{
glViewport(ctx->viewport[0], ctx->viewport[1],
ctx->viewport[2], ctx->viewport[3]);
glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
drawContext::global()->draw();
TwDraw();
glutSwapBuffers();
glutPostRedisplay();
}
void reshape(int w, int h)
{
ctx->viewport[2] = w;
ctx->viewport[3] = h;
TwWindowSize(w, h);
}
void keyboard(unsigned char key, int x, int y)
{
static bool fullScreen = false;
static int oldw = 10, oldh = 10;
if(TwEventKeyboardGLUT(key, x, y)) return;
switch(key){
case '1': GModel::current()->mesh(1); break;
case '2': GModel::current()->mesh(2); break;
case '3': GModel::current()->mesh(3); break;
case 'f':
if(fullScreen){ glutReshapeWindow(oldw, oldh); }
else{ oldw = ctx->viewport[2]; oldh = ctx->viewport[3]; glutFullScreen(); }
fullScreen = !fullScreen;
}
glutPostRedisplay();
}
void mouseMotion(int x, int y)
{
if(TwEventMouseMotionGLUT(x, y)) return;
int w = ctx->viewport[2];
int h = ctx->viewport[3];
mousePosition currPos;
currPos.set(ctx, x, y);
if(specialKey == GLUT_ACTIVE_SHIFT){
double dx = currPos.win[0] - prevPos.win[0];
double dy = currPos.win[1] - prevPos.win[1];
if(fabs(dy) > fabs(dx)) {
double fact = (4. * fabs(dy) + h) / (double)h;
ctx->s[0] *= ((dy > 0) ? fact : 1. / fact);
ctx->s[1] = ctx->s[0];
ctx->s[2] = ctx->s[0];
clickPos.recenter(ctx);
}
}
else if(specialKey == GLUT_ACTIVE_ALT){
ctx->t[0] += (currPos.wnr[0] - clickPos.wnr[0]);
ctx->t[1] += (currPos.wnr[1] - clickPos.wnr[1]);
ctx->t[2] = 0.;
}
else{
ctx->addQuaternion
((2. * prevPos.win[0] - w) / w, (h - 2. * prevPos.win[1]) / h,
(2. * currPos.win[0] - w) / w, (h - 2. * currPos.win[1]) / h);
}
prevPos.set(ctx, x, y);
glutPostRedisplay();
}
void mousePassiveMotion(int x, int y)
{
if(TwEventMouseMotionGLUT(x, y)) return;
std::vector<GVertex*> vertices;
std::vector<GEdge*> edges;
std::vector<GFace*> faces;
std::vector<GRegion*> regions;
std::vector<MElement*> elements;
bool ret = ctx->select(ENT_ALL, false, false, x, y, 5, 5,
vertices, edges, faces, regions, elements);
if(ret){
GEntity *ge = 0;
if(vertices.size()) ge = vertices[0];
else if(edges.size()) ge = edges[0];
else if(faces.size()) ge = faces[0];
else if(regions.size()) ge = regions[0];
MElement *me = elements.size() ? elements[0] : 0;
printf("%s %s\n", ge ? ge->getInfoString().c_str() : "",
me ? me->getInfoString().c_str() : "");
}
}
void mouseButton(int button, int state, int x, int y)
{
if(TwEventMouseButtonGLUT(button, state, x, y)) return;
specialKey = glutGetModifiers();
clickPos.set(ctx, x, y);
prevPos.set(ctx, x, y);
}
// AntTweakBar callbacks
void TW_CALL SetLightDirCB(const void *value, void *clientData)
{
const double *dir = (const double *)(value);
GmshSetOption("General", "Light0X", -dir[0]);
GmshSetOption("General", "Light0Y", -dir[1]);
GmshSetOption("General", "Light0Z", -dir[2]);
}
void TW_CALL GetLightDirCB(void *value, void *clientData)
{
double *dir = (double*)(value);
GmshGetOption("General", "Light0X", dir[0]);
GmshGetOption("General", "Light0Y", dir[1]);
GmshGetOption("General", "Light0Z", dir[2]);
dir[0] *= -1; dir[1] *= -1; dir[2] *= -1;
}
void TW_CALL SetInt32CB(const void *value, void *clientData)
{
int b = *(const int*)(value);
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
GmshSetOption(s.substr(0, idot), s.substr(idot + 1), (double)b);
}
void TW_CALL GetInt32CB(void *value, void *clientData)
{
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
double tmp;
GmshGetOption(s.substr(0, idot), s.substr(idot + 1), tmp);
*(int*)(value) = (int)tmp;
}
void TW_CALL SetDoubleCB(const void *value, void *clientData)
{
double b = *(const double*)(value);
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
GmshSetOption(s.substr(0, idot), s.substr(idot + 1), b);
}
void TW_CALL GetDoubleCB(void *value, void *clientData)
{
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
GmshGetOption(s.substr(0, idot), s.substr(idot + 1), *(double*)(value));
}
void TW_CALL SetColorCB(const void *value, void *clientData)
{
unsigned int b = *(const unsigned int*)(value);
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
GmshSetOption(s.substr(0, idot), s.substr(idot + 1), b);
}
void TW_CALL GetColorCB(void *value, void *clientData)
{
std::string s((const char *)clientData);
int idot = s.find_first_of('.');
GmshGetOption(s.substr(0, idot), s.substr(idot + 1), *(unsigned int*)(value));
}
void TW_CALL MenuCB(void *clientData)
{
printf("menu '%s'\n", (const char*)clientData);
}
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GmshSetOption("General", "Terminal", 1.);
GmshSetOption("View", "IntervalsType", 1.);
GmshSetOption("View", "AdaptVisualizationGrid", 1.);
GmshSetOption("View", "TargetError", 0.00001);
GmshSetOption("View", "MaxRecursionLevel", 3.);
for(int i = 1; i < argc; i++) GmshMergeFile(argv[i]);
ctx = new drawContext();
drawContext::setGlobal(new drawContextTw);
if(!TwInit(TW_OPENGL, NULL)){
printf("AntTweakBar initialization failed: %s\n", TwGetLastError());
return 1;
}
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(ctx->viewport[2], ctx->viewport[3]);
glutInitWindowPosition(100, 100);
glutCreateWindow("Gmsh Viewer");
glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutMouseFunc(mouseButton);
glutMotionFunc(mouseMotion);
glutPassiveMotionFunc(mousePassiveMotion);
glutKeyboardFunc(keyboard);
glutSpecialFunc((GLUTspecialfun)TwEventSpecialGLUT);
TwGLUTModifiersFunc(glutGetModifiers);
TwBar *bar = TwNewBar("Options");
TwDefine("Options size='200 400' color='50 50 50' alpha=128");
{
TwEnumVal axesEV[6] = { {0, "None"}, {1, "Simple axes"}, {2, "Box"},
{3, "Full grid"}, {4, "Open grid"}, {5, "Ruler"} };
TwType axesType = TwDefineEnum("AxesType", axesEV, 6);
TwAddVarCB(bar, "Axes", axesType, SetInt32CB, GetInt32CB,
(void*)"General.Axes", "group='General' help='Change axes.' ");
TwAddVarCB(bar, "LightDir", TW_TYPE_DIR3D, SetLightDirCB, GetLightDirCB,
0, "group='General' label='Light direction' close help='Change "
"light direction.' ");
{
TwAddVarCB(bar, "Background", TW_TYPE_COLOR32, SetColorCB, GetColorCB,
(void*)"General.Background", "group='GeneralColor' "
"label='Background color' ");
TwAddVarCB(bar, "BackgroundGradient", TW_TYPE_COLOR32, SetColorCB, GetColorCB,
(void*)"General.BackgroundGradient", "group='GeneralColor' "
"label='Background gradient color' ");
TwDefine("Options/GeneralColor label='Colors' close group='General' ");
}
TwDefine("Options/General close ");
}
{
TwAddVarCB(bar, "Points", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Geometry.Points", "group='Geometry' help='Draw points.' ");
TwAddVarCB(bar, "Lines", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Geometry.Lines", "group='Geometry' help='Draw lines.' ");
TwAddVarCB(bar, "Surfaces", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Geometry.Surfaces", "group='Geometry' help='Draw surfaces.' ");
TwAddVarCB(bar, "Volumes", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Geometry.Volumes", "group='Geometry' help='Draw volumes.' ");
}
{
TwAddVarCB(bar, "Vertices", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.Points", "group='Mesh' help='Draw mesh vertices.' ");
TwAddVarCB(bar, "MeshLines", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.Lines", "group='Mesh' label='Lines' help='Draw line mesh.' ");
TwAddVarCB(bar, "SurfaceEdges", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.SurfaceEdges", "group='Mesh' label='Surface edges' "
"help='Draw surface mesh edges.' ");
TwAddVarCB(bar, "SurfaceFaces", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.SurfaceFaces", "group='Mesh' label='Surface faces' "
"help='Draw surface mesh faces.' ");
TwAddVarCB(bar, "VolumeEdges", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.VolumeEdges", "group='Mesh' label='Volume edges' "
"help='Draw volume mesh edges.' ");
TwAddVarCB(bar, "VolumeFaces", TW_TYPE_BOOL32, SetInt32CB, GetInt32CB,
(void*)"Mesh.VolumeFaces", "group='Mesh' label='Volume faces' "
"help='Draw volume mesh faces.' ");
TwAddVarCB(bar, "Explode", TW_TYPE_DOUBLE, SetDoubleCB, GetDoubleCB,
(void*)"Mesh.Explode", "group='Mesh' label='Explode factor' "
"min=0 max=1 step=0.01 help='Explode mesh.' ");
TwAddVarCB(bar, "SizeFactor", TW_TYPE_DOUBLE, SetDoubleCB, GetDoubleCB,
(void*)"Mesh.CharacteristicLengthFactor", "group='Mesh' "
"label='Element size factor' min=0.01 max=100 step=0.01 ");
}
TwBar *menubar = TwNewBar("Menu");
TwDefine("Menu size='200 400' position='500 30' iconified='true' ");
TwAddButton(menubar, "Elementary entities", MenuCB, (void*)"Elementary", 0);
TwAddButton(menubar, "Physical groups", MenuCB, (void*)"Physical", 0);
TwAddButton(menubar, "Edit", MenuCB, (void*)"Edit", 0);
TwAddButton(menubar, "Reload", MenuCB, (void*)"Reload", 0);
glutMainLoop();
GmshFinalize();
return 0;
}
// lx ly lz rmax levels refcs
//
// plaqueEp.stp 0.2 0.2 0.2 0.3 3
// plaqueEpRotated.stp 0.3 0.3 0.3 0.3 3
// jonction_collee2.stp 6 6 6 10 3
// panneau_raidi_simple.stp 3 3 50 5 2 0
// plaque_trouee.stp 1 1 1 2 3
#include "Gmsh.h"
#include "GModel.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "cartesian.h"
static void insertActiveCells(double x, double y, double z, double rmax,
cartesianBox<double> &box)
{
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.getCellIJK(id1, i1, j1, k1);
int i2, j2, k2;
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));
}
static void computeLevelset(GModel *gm, cartesianBox<double> &box)
{
// tolerance for desambiguation
const double tol = box.getLC() * 1.e-12;
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 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++){
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++){
// FIXME: if there is an ambiguity assume we are inside (to
// avoid holes in the structure). This is definitely just a
// hack, as it could create pockets of matter outside the
// structure...
if(dist[j] * localdist[j] < 0 &&
fabs(fabs(dist[j]) - fabs(localdist[j])) < tol){
dist[j] = std::max(dist[j], localdist[j]);
}
else{
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());
}
static 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);
points.push_back(SPoint3(p.x(), p.y(), p.z()));
}
}
static int removeBadChildCells(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 removeBadChildCells(child);
}
static void removeParentCellsWithChildren(cartesianBox<double> *box)
{
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());
}
static 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->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;
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] [refcs=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"
" 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(" 'refcs' selects if curved surfaces should be refined\n");
return -1;
}
GmshInitialize();
GmshSetOption("General", "Terminal", 1.);
GmshMergeFile(argv[1]);
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 refineCurvedSurfaces = (argc > 7) ? atof(argv[7]) : 1;
// minimum distance between points in the cloud at the coarsest
// level
double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
// radius of the "tube" created around parts to refine at the
// coarsest level
double rtube = std::max(lx, std::max(ly, lz)) * 2.;
GModel *gm = GModel::current();
std::vector<SPoint3> points;
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> refinePoints;
if(levels > 1){
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, s, refinePoints);
// FIXME: refine this by computing e.g. "mean" curvature
if(refineCurvedSurfaces){
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 < 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;
int NZ = range.z() / lz;
if(NX < 2) NX = 2;
if(NY < 2) NY = 2;
if(NZ < 2) NZ = 2;
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());
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),
SVector3(0, 0, range.z()),
NX, NY, NZ, levels);
Msg::Info("Inserting active cells in the cartesian grid");
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())){
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;
}
// remove child cells that do not entirely fill parent cell or for
// which there is no parent neighbor; then remove parent cells that
// have children
Msg::Info("Removing cells to match X-FEM mesh topology constraints");
removeBadChildCells(&box);
removeParentCellsWithChildren(&box);
// 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);
Msg::Info("Removing cells outside the structure");
removeOutsideCells(&box);
Msg::Info("Renumbering mesh vertices across levels");
box.renumberNodes();
bool decomposeInSimplex = false;
box.writeMSH("yeah.msh", decomposeInSimplex);
Msg::Info("Done!");
GmshFinalize();
}
#include <cstring>
#include "Gmsh.h"
#include "elasticitySolver.h"
#include "PView.h"
#include "PViewData.h"
void Info (int i, char* c){
printf("%d %s\n",i,c);
}
/*
void InitializeOnelab(std::string sockName, std::string modelName) {
if (!sockName.size()) return;
loader = new onelab::remoteNetworkClient("loader", sockName);
std::vector<std::string> choices;
std::vector<onelab::string> ps;
loader->get(ps,"Elasticity/9Compute"); // ??
ps.resize(1);
ps[0].setName("Elasticity/9Compute");
ps[0].setValue("-solve -pos");
choices.push_back("-solve -pos");
choices.push_back("-pos");
choices.push_back("-solve");
ps[0].setChoices(choices);
loader->set(ps[0]);
ps.resize(1);
ps[0].setName("Elasticity/1ModelName");
ps[0].setValue(modelName);
loader->set(ps[0]);
}
void AddOnelabNumberChoice(std::string name, double val, std::string help)
{
std::vector<double> choices;
std::vector<onelab::number> ps;
loader->get(ps, name);
if(ps.size()){
choices = ps[0].getChoices();
}
else{
ps.resize(1);
ps[0].setName(name);
}
ps[0].setValue(val);
choices.push_back(val);
ps[0].setChoices(choices);
ps[0].setHelp(help);
ps[0].setShortHelp(help);
loader->set(ps[0]);
}
void AddOnelabMaterials (elasticitySolver &e) {
std::vector<onelab::number> ps;
for (int i=0;i<e.elasticFields.size();i++){
std::stringstream ss;//create a stringstream
ss << i;//add number to the stream
std::string mat = "Elasticity/LinearElasticMaterial "+ss.str()+"/Young/";
std::string help = "Young Modulus";
AddOnelabNumberChoice (mat,e.elasticFields[i]._E,help);
mat = "Elasticity/LinearElasticMaterial "+ss.str()+"/Poisson/";
help = "Poisson Ratio";
AddOnelabNumberChoice (mat,e.elasticFields[i]._nu,help);
}
}
void GetSetLoads (elasticitySolver &e) {
// todo
}
void GetSetFixations (elasticitySolver &e) {
// todo
}
*/
void WhatToDoNow(int argc, char *argv[], int &solve, std::string &fileName)
{
int i = 1;
solve = 1;
std::string sockName = "";
while (i < argc) {
if (argv[i][0] == '-') {
if (!strcmp(argv[i]+1, "solve") ||
!strcmp(argv[i]+1, "-solve")) {
solve = 1;
i++;
}
else if (!strcmp(argv[i]+1, "pos") ||
!strcmp(argv[i]+1, "-pos")) {
i++;
}
// else if (!strcmp(argv[i]+1, "onelab")) {
// i++;
// if (i<argc && argv[i][0]!='-') {
// printf("INITIALIZING SOCKET %s\n",argv[i]);
// sockName=argv[i];
// InitializeOnelab(sockName,modelName); i++;
// }
// else {
// printf("Error : Missing address of onelab server");
// }
// }
else if (!strcmp(argv[i]+1, "a")) {
i++;
exit(1);
}
else if (!strcmp(argv[i]+1, "help") || !strcmp(argv[i]+1, "h") ||
!strcmp(argv[i]+1, "-help") || !strcmp(argv[i]+1, "-h")) {
i++;
Info(0, argv[0]);
}
else if (!strcmp(argv[i]+1, "version") ||
!strcmp(argv[i]+1, "-version")) {
i++;
Info(1, argv[0]);
}
else if (!strcmp(argv[i]+1, "info") ||
!strcmp(argv[i]+1, "-info")) {
i++;
Info(2, argv[0]);
}
}
else {
std::string modelName = std::string(argv[i]);
fileName=modelName.substr(0,modelName.find_last_of(".")); // remove extension
// if (modelName)
// modelName = modelName + std::string(".fuk");
i++;
}
}
}
int main (int argc, char* argv[]){
char* a[10];
char name[245];
a[0] = name;
GmshInitialize(1, a);
GmshSetOption("General","Terminal",1.);
// for (int i=0;i<argc;i++)printf("%s\n",argv[i]);
printf("Welcome to the Simple Elasticity Program For Onelab\n");
int solve;
std::string pn;
WhatToDoNow (argc,argv, solve, pn);
elasticitySolver mySolver (1000);
mySolver.setMesh(std::string(pn+".msh").c_str());
mySolver.readInputFile(std::string(pn+".dat").c_str());
if (solve){
mySolver.solve();
PView *pvm = mySolver.buildVonMisesView("vonMises");
PView *pv = mySolver.buildDisplacementView("displacement");
// pv->getData()->writeMSH("disp.msh", false, false);
pv->getData()->writePOS("disp.pos");
pvm->getData()->writePOS("vonMises.pos");
delete pv;
}
return 0;
// solve the problem
// stop gmsh
}
//
// A simple example on how to build a GUI frontend to Gmsh using GLUT
//
#if defined(__APPLE__)
# include <GLUT/glut.h>
#else
# include <GL/glut.h>
#endif
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "drawContext.h"
drawContext *ctx = 0;
class drawContextGlut : public drawContextGlobal{
public:
void draw(){ ctx->draw3d(); ctx->draw2d(); }
const char *getFontName(int index){ return "Helvetica"; }
int getFontSize(){ return 18; }
double getStringWidth(const char *str)
{
return glutBitmapLength(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)str);
}
int getStringHeight(){ return 18; }
int getStringDescent(){ return 6; }
void drawString(const char *str)
{
for (int i = 0; i < strlen(str); i++)
glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, str[i]);
}
};
// GLUT callbacks
void display()
{
glViewport(ctx->viewport[0], ctx->viewport[1],
ctx->viewport[2], ctx->viewport[3]);
glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
drawContext::global()->draw();
glutSwapBuffers();
}
void reshape(int w, int h)
{
ctx->viewport[2] = w;
ctx->viewport[3] = h;
display();
}
void keyboard(unsigned char key, int x, int y)
{
switch(key){
case '1': GModel::current()->mesh(1); break;
case '2': GModel::current()->mesh(2); break;
case '3': GModel::current()->mesh(3); break;
}
display();
}
static int xprev = 0, yprev = 0, specialkey = 0;
void motion(int x, int y)
{
int w = ctx->viewport[2];
int h = ctx->viewport[3];
if(specialkey == GLUT_ACTIVE_SHIFT){
double dx = x - xprev;
double dy = y - yprev;
if(fabs(dy) > fabs(dx)) {
double fact = (4. * fabs(dy) + h) / (double)h;
ctx->s[0] *= ((dy > 0) ? fact : 1. / fact);
ctx->s[1] = ctx->s[0];
ctx->s[2] = ctx->s[0];
}
}
else{
ctx->addQuaternion((2. * xprev - w) / w, (h - 2. * yprev) / h,
(2. * x - w) / w, (h - 2. * y) / h);
}
xprev = x;
yprev = y;
display();
}
void mouse(int button, int state, int x, int y)
{
specialkey = glutGetModifiers();
xprev = x;
yprev = y;
}
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GmshSetOption("General", "Terminal", 1.);
GmshSetOption("View", "IntervalsType", 1.);
GmshSetOption("View", "AdaptVisualizationGrid", 1.);
GmshSetOption("View", "TargetError", 0.00001);
GmshSetOption("View", "MaxRecursionLevel", 3.);
for(int i = 1; i < argc; i++) GmshMergeFile(argv[i]);
ctx = new drawContext();
drawContext::setGlobal(new drawContextGlut);
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(ctx->viewport[2], ctx->viewport[3]);
glutInitWindowPosition(100, 100);
glutCreateWindow("GLUT Gmsh Viewer");
glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutKeyboardFunc(keyboard);
glutMotionFunc(motion);
glutMouseFunc(mouse);
glutMainLoop();
GmshFinalize();
return 0;
}
// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
//
// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
//
#include <stdio.h>
#include <sstream>
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "CellComplex.h"
#include "Homology.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GModel *m = new GModel();
m->readGEO("model.geo");
m->mesh(3);
// OR
// m->readMSH("model.msh");
// List of physical regions as domain for homology computation
// (relative to subdomain).
std::vector<int> domain;
std::vector<int> subdomain;
std::vector<int> physicalIm;
// initialize
Homology* homology = new Homology(m, domain, subdomain, physicalIm);
// find homology basis elements
homology->findHomologyBasis();
// find cohomology basis elements
homology->findCohomologyBasis();
// add 1 and 2 dimensional result chains to model
homology->addChainsToModel(1);
homology->addChainsToModel(2);
homology->addCochainsToModel(1);
homology->addCochainsToModel(2);
// write mesh with (co)homology computation result chains
m->writeMSH("model_hom.msh");
delete homology;
delete m;
GmshFinalize();
}
#include "Gmsh.h"
#include "GModel.h"
#include "MVertex.h"
#include "PView.h"
#include "PViewData.h"
#include "PluginManager.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GmshSetOption("General", "Terminal", 1.);
// load a geometry and mesh it
GModel *m = new GModel();
m->readGEO("../../tutorial/t1.geo");
m->mesh(2);
// create a node-based post-processing dataset
std::vector<GEntity*> entities;
m->getEntities(entities);
std::map<int, std::vector<double> > d;
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshVertices(); j++){
MVertex *v = entities[i]->getMeshVertex(j);
d[v->getNum()].push_back(v->x());
}
}
PView *p = new PView("f(x,y,z) = x", "NodeData", m, d);
p->getData()->writeMSH("f.msh");
// use a plugin on the dataset
PluginManager::instance()->setPluginOption("CutPlane", "A", 0.);
PluginManager::instance()->setPluginOption("CutPlane", "B", 1.);
PluginManager::instance()->setPluginOption("CutPlane", "C", 0.);
PluginManager::instance()->setPluginOption("CutPlane", "D", -0.05);
PluginManager::instance()->setPluginOption("CutPlane", "View", 0.);
PluginManager::instance()->action("CutPlane", "Run", 0);
PView::list.back()->getData()->writeMSH("fcut.msh");
GmshFinalize();
}
#include <stdio.h>
#include "Gmsh.h"
#include "GModel.h"
#include "MTetrahedron.h"
#include "MTriangle.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GModel *m = new GModel();
m->readMSH("../problem.msh");
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
GRegion *r = *it;
for(unsigned int i = 0; i < r->tetrahedra.size(); i++)
printf(" tet %d: vertices %d %d %d %d\n", r->tetrahedra[i]->getNum(),
r->tetrahedra[i]->getVertex(0)->getNum(),
r->tetrahedra[i]->getVertex(1)->getNum(),
r->tetrahedra[i]->getVertex(2)->getNum(),
r->tetrahedra[i]->getVertex(3)->getNum());
}
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
GFace *f = *it;
for(unsigned int i = 0; i < f->triangles.size(); i++)
printf(" tri %d: vertices %d %d %d\n", f->triangles[i]->getNum(),
f->triangles[i]->getVertex(0)->getNum(),
f->triangles[i]->getVertex(1)->getNum(),
f->triangles[i]->getVertex(2)->getNum());
}
delete m;
GmshFinalize();
}
#include <stdio.h>
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "discreteRegion.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GmshSetOption("General", "Terminal", 1.);
GmshSetOption("General", "Verbosity", 4.);
GmshSetOption("Mesh", "CharacteristicLengthExtendFromBoundary", 1.);
GmshSetOption("Mesh", "OldRefinement", 1.);
GmshSetOption("Mesh", "CharacteristicLengthMin", 0.1);
GmshSetOption("Mesh", "CharacteristicLengthMax", 0.1);
GmshSetOption("Mesh", "Optimize", 0.); // not yet: need boundary!
GModel *m = new GModel();
m->readMSH("cube.msh");
// discreteRegion *gr = new discreteRegion(m);
// MVertex *v0 = new MVertex(x, y, z, gr, tag);
// MVertex *v1 = new MVertex(x, y, z, gr, tag);
// ...
// gr->mesh_vertices.push_back(v0);
// ...
// MTetrahedron *t = new MTetrahedron(v0, v1, v2, v3, tag);
// gr->tetrahedra.push_back(t);
//MVertex *v2 = gr->mesh_vertices[2];
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
discreteRegion *r = dynamic_cast<discreteRegion*>(*it);
if(r){
printf("found a discrete region (%d) with %d tetrahedra\n", r->tag(), r->tetrahedra.size());
r->remesh();
printf("after remeshing: %d tetrahedra\n", r->tetrahedra.size());
}
}
m->writeMSH("new.msh");
delete m;
GmshFinalize();
}
#include <stdio.h>
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "VertexArray.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
GmshSetOption("Mesh", "Algorithm", 5.);
GmshSetOption("General", "Terminal", 1.);
GModel *m = new GModel();
m->readMSH("bunny.msh");
m->fillVertexArrays();
std::vector<GEntity*> entities;
m->getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
printf("coucou entite %d (dimension %d)\n", ge->tag(), ge->dim());
if(ge->va_triangles)
printf(" j'ai un va de triangles: %d vertex\n", ge->va_triangles->getNumVertices());
if(ge->va_lines)
printf(" j'ai un va de lignes: %d vertex\n", ge->va_lines->getNumVertices());
}
delete m;
GmshFinalize();
}
onelab.number Material1.Create(1);
onelab.string Material.Add(Material1); Material.Add(Material2); Material.Add(Material3);
onelab.iftrue(Material1)
onelab.number E.Create(100.e9); Nu.Create(0.3);
onelab.endif
onelab.iftrue(Material2)
onelab.number E.Create(200.e9); Nu.Create(0.25);
onelab.endif
onelab.iftrue(Material3)
onelab.number E.Create(80.e9); Nu.Create(0.5);
onelab.endif
ElasticDomain 7 onelab.getValue(E) onelab.getValue(Nu);
EdgeDisplacement 8 0 0
EdgeDisplacement 8 1 0
EdgeDisplacement 8 2 0
EdgeForce 9 100e6 0 0
Mesh.RecombinationAlgorithm=1;
Mesh.Algorithm = 8;
unit = 1.0e-02 ;
DefineConstant[ H = {4.5 * unit, Min 1 *unit, Max 8.5 *unit, Step 1*unit,
Name "Parameters/Geometry/Beam Height"} ] ;
DefineConstant[ L = {20 * unit, Min 10 *unit, Max 200 *unit, Step 1*unit,
Name "Parameters/Geometry/Beam Width"} ] ;
DefineConstant[ lc = {2 * unit, Min .1 *unit, Max 10 *unit, Step .1*unit,
Name "Parameters/Geometry/Mesh Size"} ] ;
Point(1) = {0, 0, 0, lc};
Point(2) = {H, 0, 0, lc};
Point(3) = {H, L, 0, lc};
Point(4) = {0, L, 0, lc};
Line(1) = {4, 3};
Line(2) = {3, 2};
Line(3) = {2, 1};
Line(4) = {1, 4};
Line Loop(5) = {2, 3, 4, 1};
Plane Surface(6) = {5};
Physical Surface(7) = {6};
Physical Line(8) = {1};
Physical Line(9) = {3};
Recombine Surface {6};
File moved
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment