diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6c6831fde44ed83fc2d3773102a68812dcf5cdf3..9f1620c42c4767d1f9236568eb0062fad2c2b448 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1330,7 +1330,7 @@ bool GFace::fillPointCloud(double maxDist,
SPoint2 p = p0 * (1. - u - v) + p1 * u + p2 * v;
GPoint gp(point(p));
points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
- if (uvpoints)uvpoints->push_back(p);
+ if(uvpoints) uvpoints->push_back(p);
if(normals) normals->push_back(normal(p));
}
}
@@ -1348,8 +1348,8 @@ bool GFace::fillPointCloud(double maxDist,
double t2 = b2.low() + v * (b2.high() - b2.low());
GPoint gp = point(t1, t2);
points->push_back(SPoint3(gp.x(), gp.y(), gp.z()));
- if (uvpoints)uvpoints->push_back(SPoint2(t1,t2));
- if(normals) normals->push_back(normal(SPoint2(t1,t2)));
+ if(uvpoints) uvpoints->push_back(SPoint2(t1, t2));
+ if(normals) normals->push_back(normal(SPoint2(t1, t2)));
}
}
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index afc87653000a2322bedfb6dc0ae4784bf9a0d52e..938f29a521adba5d67b5a5fd6f1bd597f98c3126 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -263,7 +263,7 @@ class GFace : public GEntity
// points are at most maxDist apart
bool fillPointCloud(double maxDist,
std::vector<SPoint3> *points,
- std::vector<SPoint2> *uvpoints,
+ std::vector<SPoint2> *uvpoints=0,
std::vector<SVector3> *normals=0);
// apply Lloyd's algorithm to the mesh
diff --git a/utils/api_demos/CMakeLists.txt b/utils/api_demos/CMakeLists.txt
index 97db72d7ce346335df436bb87a86eb1c25a8d54a..b10031326e7f059ee6b9469c322f898e7cb7e4cb 100644
--- a/utils/api_demos/CMakeLists.txt
+++ b/utils/api_demos/CMakeLists.txt
@@ -14,8 +14,10 @@ project(api_demos CXX)
add_subdirectory(../.. "${CMAKE_CURRENT_BINARY_DIR}/gmsh")
include_directories(../../Common ../../Numeric ../../Geo ../../Mesh
- ../../Solver ../../Post ../../Plugin ../../Graphics ../../contrib/ANN/include
- ../../contrib/DiscreteIntegration ${GMSH_EXTERNAL_INCLUDE_DIRS}
+ ../../Solver ../../Post ../../Plugin ../../Graphics
+ ../../contrib/ANN/include ../../contrib/MathEx ../../contrib/kbipack
+ ../../contrib/DiscreteIntegration
+ ${GMSH_EXTERNAL_INCLUDE_DIRS}
${CMAKE_CURRENT_BINARY_DIR}/gmsh/Common)
if(APPLE)
@@ -37,14 +39,11 @@ add_executable(mainElasticity mainElasticity.cpp)
target_link_libraries(mainElasticity shared)
add_executable(mainGlut mainGlut.cpp)
-target_link_libraries(mainGlut lib ${GMSH_EXTERNAL_LIBRARIES} ${glut})
+target_link_libraries(mainGlut shared ${glut})
add_executable(mainHomology mainHomology.cpp)
target_link_libraries(mainHomology shared)
-add_executable(mainLevelset mainLevelset.cpp)
-target_link_libraries(mainLevelset shared)
-
add_executable(mainOcc mainOcc.cpp)
target_link_libraries(mainOcc shared)
diff --git a/utils/api_demos/mainGlut.cpp b/utils/api_demos/mainGlut.cpp
index 1d17c40538e0bebc2ccc36246e5608771fde4aec..67cb5f177dfcb1b25f8a62af4afbb129cf4eb35d 100644
--- a/utils/api_demos/mainGlut.cpp
+++ b/utils/api_demos/mainGlut.cpp
@@ -1,4 +1,3 @@
-
//
// A simple example on how to build a GUI frontend to Gmsh using GLUT
//
@@ -6,133 +5,43 @@
#if defined(__APPLE__)
# include <GLUT/glut.h>
#else
-# include <GL/gl.h>
# include <GL/glut.h>
-# include <GL/glu.h>
#endif
#include "Gmsh.h"
-#include <stdio.h>
-#include "string.h"
#include "GModel.h"
#include "MElement.h"
-#include "Context.h"
#include "drawContext.h"
-#include "Trackball.h"
-#include "Camera.h"
-
-
-typedef struct {
- double r,g,b;
-} COLOUR;
-typedef struct {
- unsigned char r,g,b,a;
-} PIXELA;
-using namespace std;
drawContext *ctx = 0;
-Camera *camera;
-mouseAndKeyboard mouseandkeys;
-
-
-static int xprev = 0, yprev = 0, specialkey = 0;
class drawContextGlut : public drawContextGlobal{
-public:
+ public:
void draw(){ ctx->draw3d(); ctx->draw2d(); }
const char *getFontName(int index){ return "Helvetica"; }
- int getFontSize(){ return 12; }
+ int getFontSize(){ return 18; }
double getStringWidth(const char *str)
{
- return glutBitmapLength(GLUT_BITMAP_HELVETICA_12, (const unsigned char*)str);
+ return glutBitmapLength(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)str);
}
- int getStringHeight(){ return 12; }
+ 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]);
+ for (int i = 0; i < strlen(str); i++)
+ glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, str[i]);
}
-
};
-
-
// GLUT callbacks
void display()
{
-
- if (!camera->stereoEnable) {
- glViewport(ctx->viewport[0], ctx->viewport[1],
- ctx->viewport[2], ctx->viewport[3]);
- glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
- camera->giveViewportDimension(ctx->viewport[2],ctx->viewport[3]);
-
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- glFrustum(camera->glFleft,camera->glFright,camera->glFbottom,camera->glFtop,camera->glFnear,camera->glFfar);
- glMatrixMode(GL_MODELVIEW);
- glDrawBuffer(GL_BACK);
- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
- glLoadIdentity();
- gluLookAt(camera->position.x, camera->position.y, camera->position.z,
- camera->target.x, camera->target.y, camera->target.z,
- camera->up.x, camera->up.y, camera->up.z);
- ctx->draw3d();
- ctx->draw2d();
- glutSwapBuffers();
-
- }
- else {
- // cout<<" display3D();"<<endl;
- glViewport(ctx->viewport[0], ctx->viewport[1],
- ctx->viewport[2], ctx->viewport[3]);
- glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
- camera->giveViewportDimension(ctx->viewport[2],ctx->viewport[3]);
-
- //right eye
- XYZ eye =camera->eyesep / 2.0* camera->right;
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- double left = - camera->screenratio * camera->wd2 - 0.5 * camera->eyesep * camera->ndfl;
- double right = camera->screenratio * camera->wd2 - 0.5 * camera->eyesep * camera->ndfl;
- double top = camera->wd2;
- double bottom = - camera->wd2;
- glFrustum(left,right,bottom,top,camera->glFnear,camera->glFfar);
- glMatrixMode(GL_MODELVIEW);
- glDrawBuffer(GL_BACK_RIGHT);
- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
- glLoadIdentity();
- gluLookAt(camera->position.x+eye.x, camera->position.y+eye.y, camera->position.z+eye.z,
- camera->target.x+eye.x, camera->target.y+eye.y, camera->target.z+eye.z,
- camera->up.x, camera->up.y, camera->up.z);
- ctx->draw3d();
- ctx->draw2d();
-
- //left eye
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- left = - camera->screenratio * camera->wd2 + 0.5 * camera->eyesep * camera->ndfl;
- right = camera->screenratio * camera->wd2 + 0.5 * camera->eyesep * camera->ndfl;
- top = camera->wd2;
- bottom = - camera->wd2;
- glFrustum(left,right,bottom,top,camera->glFnear,camera->glFfar);
-
- glMatrixMode(GL_MODELVIEW);
- glDrawBuffer(GL_BACK_LEFT);
- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
- glLoadIdentity();
- gluLookAt(camera->position.x-eye.x, camera->position.y-eye.y, camera->position.z-eye.z,
- camera->target.x-eye.x, camera->target.y-eye.y, camera->target.z-eye.z,
- camera->up.x, camera->up.y, camera->up.z);
- ctx->draw3d();
- ctx->draw2d();
-
- glutSwapBuffers();
- }
-
+ 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;
@@ -140,239 +49,74 @@ void reshape(int w, int 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();
+}
-void keyboard(unsigned char key, int x, int y) { display();}
-
+static int xprev = 0, yprev = 0, specialkey = 0;
void motion(int x, int y)
{
- int w = (ctx->viewport[2] - ctx->viewport[0]);
- int h = (ctx->viewport[3] - ctx->viewport[1]);
- //rotate
- if (mouseandkeys.button_left_down && (mouseandkeys.mode!= GLUT_ACTIVE_CTRL) ){
- double x_r = 2.*(1.*x - w/2.)/w;
- double y_r = 2.*(-1.*y + h/2.)/h;
- double xprev_r=2.*(1.*xprev - w/2.)/w;
- double yprev_r=2.*(-1.*yprev + h/2.)/h;
- double q[4];
- trackball(q,xprev_r,yprev_r,x_r,y_r);
- camera->rotate(q);
- xprev = x;
- yprev = y;
- }
- //zoom
- if (mouseandkeys.button_middle_down ){
- double dy= y-yprev;
- double factor =( CTX::instance()->zoomFactor * fabs(dy) +(double) h) / (double)h;
- factor = ((dy > 0) ? factor : 1./factor);
- camera->distance=fabs(1./factor*camera->ref_distance);
- camera->position.x=camera->target.x-camera->distance*camera->view.x;
- camera->position.y=camera->target.y-camera->distance*camera->view.y;
- camera->position.z=camera->target.z-camera->distance*camera->view.z;
- }
- // translate
- if (mouseandkeys.button_right_down ){
- double x_r = 2.*(1.*x - w/2.)/w;
- double y_r = 2.*(1.*y - h/2.)/h;
- double xprev_r=2.*(1.*xprev - w/2.)/w;
- double yprev_r=2.*(1.*yprev - h/2.)/h;
- double theta_x=camera->aperture*(x_r-xprev_r)*0.0174532925/2. ;
- double theta_y=camera->aperture*(y_r-yprev_r)*0.0174532925/2. ;
- camera->moveRight(theta_x);
- camera->moveUp(theta_y);
- xprev = x;
- yprev = 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];
+ }
}
- if (!mouseandkeys.button_middle_down ){
- xprev = x;
- yprev = y;
+ 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;
- if (button == GLUT_LEFT_BUTTON) {
- mouseandkeys.button_left_down=!state;
- }
- else if (button == GLUT_MIDDLE_BUTTON) {
- mouseandkeys.button_middle_down=!state;
- }
- else {
- mouseandkeys.button_right_down=!state;
- }
- camera->ref_distance=camera->distance;
-}
-
-
-
-void processSpecialKeys(int key, int x, int y) {
-
- mouseandkeys.mode = glutGetModifiers();
- if (mouseandkeys.mode == GLUT_ACTIVE_CTRL) {
- switch(key){
- case 100 : /* 'left' */
- display();
- camera->focallength=camera->focallength*.99;
- camera->eyesep=(camera->focallength/camera->distance)*camera->distance*camera->screenratio;
- break;
- case 102 : /* 'right' */
- camera->focallength=camera->focallength*1.01;
- camera->eyesep=(camera->focallength/camera->distance)*camera->distance*camera->screenratio;
- display();
- break;
- }
- }
- else{
- switch(key){
- mouseandkeys.key= key ;
- case 101 : /* 'up' */
- camera->focallength*=1.1;
- break;
- case 103 : /* 'down' */
- camera->focallength*=0.9;
- break;
- case 100 : /* 'left' */
- camera->eyesep*=.9;
- break;
- case 102 : /* 'right' */
- camera->eyesep*=1.1;
- break;
- }
- }
- camera->update();
- display();
-
- cout<<"eyesep = "<< camera->eyesep<<" / focallength = "<< camera->focallength<<" / screenratio = ";
- cout<< camera->screenratio<<" / distance = "<< camera->distance<<" / closeness = "<< camera->closeness<<endl;
- mouseandkeys.mode = 0;
-}
-
-
-
-void processNormalKeys(unsigned char key, int x, int y) {
- double def;
-
- if (key != 0) {
- mouseandkeys.key= key ;
- mouseandkeys.mode = glutGetModifiers();
- if (mouseandkeys.mode == GLUT_ACTIVE_CTRL && key==17) exit(0);
- }
- switch(key){
- case 49: GModel::current()->mesh(1); break; /* '1' */
- case 50: GModel::current()->mesh(2); break; /* '2' */
- case 51: GModel::current()->mesh(3); break; /* '3' */
- case 48 : /* '0' */
- camera->lookAtCg();
- display();
- break;
- case 114 : /* 'R' */
- camera->init();
- display();
- break;
- case 'f' : /* 'F' */
- camera->screenwidth=ctx->viewport[2];
- camera->screenheight=ctx->viewport[3];
- glutFullScreen();
- display();
- break;
- case 27 : /* 'ech' */
- glutReshapeWindow(camera->screenwidth ,camera->screenheight );
- //glutReshapeWindow(500,500 );
- display();
- break;
- case 100 : ; /* 'D' */
- GmshSetOption("View","VectorType",5.);
- display();
- break;
- case 43 : /* '+' */
- GmshGetOption("View", "DisplacementFactor", def);
- def*=2.;
- GmshSetOption("View","DisplacementFactor",def);
- display();
- break;
- case 45 : ; /* '-' */
- GmshGetOption("View", "DisplacementFactor", def);
- def*=.5;
- GmshSetOption("View","DisplacementFactor",def);
- display();
- break;
- }
-
}
-
-//*******************************************************
-//*******************************************************
-//*******************************************************
-
int main(int argc, char **argv)
{
-
GmshInitialize(argc, argv);
GmshSetOption("General", "Terminal", 1.);
- GmshSetOption("General", "Stereo", 0.);
- GmshSetOption("General", "Camera", 1.);
- GmshSetOption("General", "Orthographic", 0.);
- GmshSetOption("General", "TrackballHyperbolicSheet", 0.);
-
- if (strstr(argv[1],"-s") != NULL){
- camera->stereoEnable = true;
- cout<<"mode STEREO"<<endl;
- GmshSetOption("General", "Stereo", 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);
- camera= &(ctx->camera);
- camera->init();
glutInit(&argc, argv);
- if (camera->stereoEnable) {
- glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH | GLUT_STEREO);
- }
- else {
- glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
- }
+ glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(ctx->viewport[2], ctx->viewport[3]);
- glutInitWindowPosition(400,10);
- glutInitWindowSize(800,800);
+ glutInitWindowPosition(100, 100);
glutCreateWindow("GLUT Gmsh Viewer");
glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutKeyboardFunc(keyboard);
glutMotionFunc(motion);
glutMouseFunc(mouse);
- glutKeyboardFunc(processNormalKeys);
- glutSpecialFunc(processSpecialKeys);
-
- cout<<"-------------------------------------"<<endl;
- cout<<" SHORTCUTS "<<endl;
- cout<<"-------------------------------------"<<endl;
- cout<<" Ctrl+Q Quit "<<endl;
- cout<<" 1 mesh line "<<endl;
- cout<<" 2 mesh surface "<<endl;
- cout<<" 3 mesh volume "<<endl;
- cout<<" R resize "<<endl;
- cout<<" 0 origine "<<endl;
- cout<<" option -s stereo on "<<endl;
- cout<<" F full screen on "<<endl;
- cout<<" ech full screen off"<<endl;
- cout<<" D displacement field"<<endl;
- cout<<" + deformed +"<<endl;
- cout<<" - deformed -"<<endl;
- cout<<" up focal length +"<<endl;
- cout<<" down focal length -"<<endl;
- cout<<" left eye sep -"<<endl;
- cout<<" right eye sep +"<<endl;
- cout<<"-------------------------------------"<<endl;
glutMainLoop();
+
GmshFinalize();
return 0;
}
diff --git a/utils/api_demos/mainHomology.cpp b/utils/api_demos/mainHomology.cpp
index c6b324e27fd432b3a330077a208b584686f9d7f9..1839ac8070b431252fb3ac86d10f9a7920b4787f 100644
--- a/utils/api_demos/mainHomology.cpp
+++ b/utils/api_demos/mainHomology.cpp
@@ -28,9 +28,10 @@ int main(int argc, char **argv)
// (relative to subdomain).
std::vector<int> domain;
std::vector<int> subdomain;
+ std::vector<int> physicalIm;
// initialize
- Homology* homology = new Homology(m, domain, subdomain);
+ Homology* homology = new Homology(m, domain, subdomain, physicalIm);
// find homology basis elements
homology->findHomologyBasis();
diff --git a/utils/api_demos/mainLevelset.cpp b/utils/api_demos/mainLevelset.cpp
deleted file mode 100644
index d6ca525fa33cf2aa050f6783dffccf96c0a71198..0000000000000000000000000000000000000000
--- a/utils/api_demos/mainLevelset.cpp
+++ /dev/null
@@ -1,701 +0,0 @@
-
-#include <time.h>
-#include <iostream>
-#include <queue>
-#include <limits>
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "MHexahedron.h"
-#include "MTetrahedron.h"
-#include "MQuadrangle.h"
-#include "MTriangle.h"
-#include "MLine.h"
-#include "Integration3D.h"
-#include "DILevelset.h"
-
-#define PI 3.14159265
-
-// ----------------------------------------------------------------------
-// --------------------------OUTPUT FILE --------------------------------
-// ----------------------------------------------------------------------
-
-int nbZeros(double number){
- if (number != 0) {
- double lg = log10(fabs(number));
- return (int)ceil(-lg) -1;
- }
- else
- return -1;
-}
-
-int digits(double number){
- int nC=0; double n=number, intpart; //printf("%-#15.15g \n",n);
- double eps = 1.e-15; //printf("%-#15.15g->%-#15.15g>%-#15.15g nC=%d\n",n,modf(n,&intpart),eps,nC);
- while(fabs(modf(n,&intpart))-1. < -eps && fabs(modf(n,&intpart)) > eps){
- nC++; n *=10; eps *=10; //printf("%-#15.17g->%-#15.17g>%-#15.17g nC=%d\n",n,modf(n,&intpart),eps,nC);
- } //printf("\n\n");
- return nC;
-}
-
-
-void writePosFile (const std::vector<DI_IntegrationPoint> &ip, const std::vector<DI_IntegrationPoint> &ipS, const std::vector<DI_CuttingPoint> &cp,
- std::vector<DI_Line> surfLines, std::vector<DI_Triangle> &surfTriangles, std::vector<DI_Quad> &surfQuads,
- const std::vector<DI_Tetra> &Tetras, const std::vector<DI_Hexa> &Hexas)
-{
- // check if a line has been created twice.
- printf("checking %d lines, %d triangles and %d quadrangles\n",surfLines.size(),surfTriangles.size(),surfQuads.size());
- for(int i=0;i<(int)surfLines.size();i++){
- for(int j=i+1;j<(int)surfLines.size();j++){
- int b=0;
- for(int k=0;k<2;k++){
- for(int l=0;l<2;l++){
- if(surfLines[i].pt(k).equal(surfLines[j].pt(l))) b++;
- }
- }
- if(b==2) {printf("BIIIP : line created twice (%d,%d)\n",i,j); surfLines[i].print(); surfLines[j].print();
- }
- }
- }
- // check if a triangle has been created twice.
- double tot=0.;
- for(int i=0;i<(int)surfTriangles.size();i++){
- tot += surfTriangles[i].surface();
- for(int j=i+1;j<(int)surfTriangles.size();j++){
- int b=0;
- for(int k=0;k<3;k++){
- for(int l=0;l<3;l++){
- if(surfTriangles[i].pt(k).equal(surfTriangles[j].pt(l))) b++;
- }
- }
- if(b==3) {printf("BIIIP : triangle created twice (%d,%d)\n",i,j); surfTriangles[i].print(); surfTriangles[j].print();
- double lsi[3] = {10,10,10};
- surfTriangles[i].addLs(lsi);
- }
- }
- }
- // check if a quad has been created twice.
- for(int i=0;i<(int)surfQuads.size();i++){
- tot +=surfQuads[i].surface();
- for(int j=i+1;j<(int)surfQuads.size();j++){
- int b=0;
- for(int k=0;k<4;k++){
- for(int l=0;l<4;l++){
- if(surfQuads[i].pt(k).equal(surfQuads[j].pt(l))) b++;
- }
- }
- if(b==4) {printf("BIIIP : quadrangle created twice (%d,%d)\n",i,j); surfQuads[i].print(); surfQuads[j].print();
- double lsi[3] = {10,10,10};
- surfQuads[i].addLs(lsi);
- }
- }
- }
- printf("somme des \"surface()\":%g\n",tot);
-
- //check if a cuttingPoint has been created twice
- for(int i=0;i<(int)cp.size();i++) {
- for(int j=i+1;j<(int)cp.size();j++){
- if(cp[i].equal(cp[j]))
- {printf("BIIIP : CuttingPoint created twice (%d,%d)\n",i,j); cp[i].print(); cp[j].print();}
- }
- }//*/
-
- // check quality of the triangles and tetras
- double worstTr=1, qMTr=0, Smin=10;
- for(int i=0;i<(int)surfTriangles.size();i++){
- double qt=surfTriangles[i].quality();
- if(qt<worstTr) worstTr=qt;
- if(surfTriangles[i].surface()<Smin) Smin=surfTriangles[i].surface();
- qMTr += qt;
- }
- double worstT=1, qMT=0;
- for(int i=0;i<(int)Tetras.size();i++){
- double qt=Tetras[i].quality();
- if(qt<worstT) worstT=qt;
- qMT += qt;
- }
- printf("Quality :\n");
- if(surfTriangles.size()) printf(" triangles : worst=%g, mean=%g\n",worstTr,qMTr/surfTriangles.size());
- if(Tetras.size()) printf(" tetras : worst=%g, mean=%g\n",worstT,qMT/Tetras.size());
- printf("Smin = %g\n",Smin);
-
- printf("Writing file\n");
- FILE *f = fopen("MeshCut.pos","w");
- fprintf(f,"View\"\"{\n");
- int nbLn = surfLines.size();
- for (int i=0;i<nbLn;i++) {
- //if(surfLines[i].lsTag()!=-1){
- /*if(surfLines[i].isQuad()){
- fprintf(f,"SL2(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- surfLines[i].x(0), surfLines[i].y(0), surfLines[i].z(0), surfLines[i].x(1), surfLines[i].y(1), surfLines[i].z(1),
- surfLines[i].x2(0), surfLines[i].y2(0), surfLines[i].z2(0),
- surfLines[i].pt(0).ls(),surfLines[i].pt(1).ls(),surfLines[i].mid(0).ls());
- }
- else{*/
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g){%d,%d};\n",//{%g,%g};\n",
- surfLines[i].x(0), surfLines[i].y(0), surfLines[i].z(0), surfLines[i].x(1), surfLines[i].y(1), surfLines[i].z(1),
- surfLines[i].lsTag(), surfLines[i].lsTag());
- //surfLines[i].pt(0).ls(),surfLines[i].pt(1).ls());
- //}//}
- }
- int nbTr = surfTriangles.size();
- for (int i=0;i<nbTr;i++) {
- //if(surfTriangles[i].lsTag()!=-1){
- //if(surfTriangles[i].x(0)>=xm && surfTriangles[i].x(0)<=xM && surfTriangles[i].y(0)>=ym && surfTriangles[i].y(0)<=yM && surfTriangles[i].z(0)>=zm && surfTriangles[i].z(0)<=zM && surfTriangles[i].x(1)>=xm && surfTriangles[i].x(1)<=xM && surfTriangles[i].y(1)>=ym && surfTriangles[i].y(1)<=yM && surfTriangles[i].z(1)>=zm && surfTriangles[i].z(1)<=zM && surfTriangles[i].x(2)>=xm && surfTriangles[i].x(2)<=xM && surfTriangles[i].y(2)>=ym && surfTriangles[i].y(2)<=yM && surfTriangles[i].z(2)>=zm && surfTriangles[i].z(2)<=zM)
- /*if(surfTriangles[i].isQuad()){
- fprintf(f,"ST2(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d,%d,%d};\n", //{%g,%g,%g,%g,%g,%g};\n",
- surfTriangles[i].x(0), surfTriangles[i].y(0), surfTriangles[i].z(0), surfTriangles[i].x(1), surfTriangles[i].y(1), surfTriangles[i].z(1),
- surfTriangles[i].x(2), surfTriangles[i].y(2), surfTriangles[i].z(2), surfTriangles[i].x2(0), surfTriangles[i].y2(0), surfTriangles[i].z2(0),
- surfTriangles[i].x2(1), surfTriangles[i].y2(1), surfTriangles[i].z2(1), surfTriangles[i].x2(2), surfTriangles[i].y2(2), surfTriangles[i].z2(2),
- surfTriangles[i].lsTag(), surfTriangles[i].lsTag(),surfTriangles[i].lsTag(),surfTriangles[i].lsTag(), surfTriangles[i].lsTag(),surfTriangles[i].lsTag());
- //surfTriangles[i].pt(0).ls(),surfTriangles[i].pt(1).ls(),surfTriangles[i].pt(2).ls(),surfTriangles[i].mid(0).ls(),surfTriangles[i].mid(1).ls(),surfTriangles[i].mid(2).ls());
- }
- else{*/
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",//{%g,%g,%g};\n",
- surfTriangles[i].x(0), surfTriangles[i].y(0), surfTriangles[i].z(0), surfTriangles[i].x(1), surfTriangles[i].y(1), surfTriangles[i].z(1),
- surfTriangles[i].x(2), surfTriangles[i].y(2), surfTriangles[i].z(2),
- //surfTriangles[i].lsTag(), surfTriangles[i].lsTag(),surfTriangles[i].lsTag());
- surfTriangles[i].pt(0).ls(),surfTriangles[i].pt(1).ls(),surfTriangles[i].pt(2).ls());
- //}//}
- }
- int nbQ = surfQuads.size();
- for (int i=0;i<nbQ;i++) {
- //if(surfQuads[i].lsTag()!=-1){
- /*if(surfQuads[i].isQuad()){
- fprintf(f,"SQ2(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
- surfQuads[i].x(0), surfQuads[i].y(0), surfQuads[i].z(0), surfQuads[i].x(1), surfQuads[i].y(1), surfQuads[i].z(1),
- surfQuads[i].x(2), surfQuads[i].y(2), surfQuads[i].z(2), surfQuads[i].x(3), surfQuads[i].y(3), surfQuads[i].z(3),
- surfQuads[i].x2(0), surfQuads[i].y2(0), surfQuads[i].z2(0), surfQuads[i].x2(1), surfQuads[i].y2(1), surfQuads[i].z2(1),
- surfQuads[i].x2(2), surfQuads[i].y2(2), surfQuads[i].z2(2), surfQuads[i].x2(3), surfQuads[i].y2(3), surfQuads[i].z2(3),
- surfQuads[i].x2(4), surfQuads[i].y2(4), surfQuads[i].z2(4),
- surfQuads[i].pt(0).ls(),surfQuads[i].pt(1).ls(),surfQuads[i].pt(2).ls(),surfQuads[i].pt(3).ls(),
- surfQuads[i].mid(0).ls(),surfQuads[i].mid(1).ls(),surfQuads[i].mid(2).ls(),surfQuads[i].mid(3).ls(),surfQuads[i].mid(4).ls());
- }
- else {*/
- fprintf(f,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n", //{%g,%g,%g,%g};\n",
- surfQuads[i].x(0), surfQuads[i].y(0), surfQuads[i].z(0), surfQuads[i].x(1), surfQuads[i].y(1), surfQuads[i].z(1),
- surfQuads[i].x(2), surfQuads[i].y(2), surfQuads[i].z(2), surfQuads[i].x(3), surfQuads[i].y(3), surfQuads[i].z(3),
- //surfQuads[i].lsTag(), surfQuads[i].lsTag(),surfQuads[i].lsTag(),surfQuads[i].lsTag());
- surfQuads[i].pt(0).ls(),surfQuads[i].pt(1).ls(),surfQuads[i].pt(2).ls(),surfQuads[i].pt(3).ls());
- //}//}
- }
- int nbT = Tetras.size();
- //double xm=3.33, xM=3.89, ym=2.22, yM=2.78, zm=1.66, zM=2.23;
- for (int i=0;i<nbT;i++) {
- //if(Tetras[i].sign()<=0){
- //if(Tetras[i].x(0)>=xm && Tetras[i].x(0)<=xM && Tetras[i].y(0)>=ym && Tetras[i].y(0)<=yM && Tetras[i].z(0)>=zm && Tetras[i].z(0)<=zM && Tetras[i].x(1)>=xm && Tetras[i].x(1)<=xM && Tetras[i].y(1)>=ym && Tetras[i].y(1)<=yM && Tetras[i].z(1)>=zm && Tetras[i].z(1)<=zM && Tetras[i].x(2)>=xm && Tetras[i].x(2)<=xM && Tetras[i].y(2)>=ym && Tetras[i].y(2)<=yM && Tetras[i].z(2)>=zm && Tetras[i].z(2)<=zM && Tetras[i].x(3)>=xm && Tetras[i].x(3)<=xM && Tetras[i].y(3)>=ym && Tetras[i].y(3)<=yM && Tetras[i].z(3)>=zm && Tetras[i].z(3)<=zM){
- /*if(Tetras[i].isQuad()){
- fprintf(f,"SS2(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
- Tetras[i].x(0), Tetras[i].y(0), Tetras[i].z(0), Tetras[i].x(1), Tetras[i].y(1), Tetras[i].z(1),
- Tetras[i].x(2), Tetras[i].y(2), Tetras[i].z(2), Tetras[i].x(3), Tetras[i].y(3), Tetras[i].z(3),
- Tetras[i].x2(0), Tetras[i].y2(0), Tetras[i].z2(0), Tetras[i].x2(3), Tetras[i].y2(3), Tetras[i].z2(3),
- Tetras[i].x2(1), Tetras[i].y2(1), Tetras[i].z2(1), Tetras[i].x2(2), Tetras[i].y2(2), Tetras[i].z2(2),
- Tetras[i].x2(4), Tetras[i].y2(4), Tetras[i].z2(4), Tetras[i].x2(5), Tetras[i].y2(5), Tetras[i].z2(5),
- Tetras[i].pt(0).ls(),Tetras[i].pt(1).ls(),Tetras[i].pt(2).ls(),Tetras[i].pt(3).ls(),
- Tetras[i].mid(0).ls(),Tetras[i].mid(1).ls(),Tetras[i].mid(2).ls(),Tetras[i].mid(3).ls(),Tetras[i].mid(4).ls(),Tetras[i].mid(5).ls());
- }
- else{*/
- fprintf(f,"SS(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",//{%g,%g,%g,%g};\n",
- Tetras[i].x(0), Tetras[i].y(0), Tetras[i].z(0), Tetras[i].x(1), Tetras[i].y(1), Tetras[i].z(1),
- Tetras[i].x(2), Tetras[i].y(2), Tetras[i].z(2), Tetras[i].x(3), Tetras[i].y(3), Tetras[i].z(3),//1.,1.,1.,1.);
- Tetras[i].lsTag(),Tetras[i].lsTag(),Tetras[i].lsTag(),Tetras[i].lsTag());
- //Tetras[i].pt(0).ls(),Tetras[i].pt(1).ls(),Tetras[i].pt(2).ls(),Tetras[i].pt(3).ls());
- //}//}
- }
- int nbH = Hexas.size();
- for (int i=0;i<nbH;i++) {
- //if(Hexas[i].sign()<=0){
- fprintf(f,"SH(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d,%d,%d,%d,%d};\n",//{%g,%g,%g,%g,%g,%g,%g,%g};\n",
- Hexas[i].x(0), Hexas[i].y(0), Hexas[i].z(0), Hexas[i].x(1), Hexas[i].y(1), Hexas[i].z(1),
- Hexas[i].x(2), Hexas[i].y(2), Hexas[i].z(2), Hexas[i].x(3), Hexas[i].y(3), Hexas[i].z(3),
- Hexas[i].x(4), Hexas[i].y(4), Hexas[i].z(4), Hexas[i].x(5), Hexas[i].y(5), Hexas[i].z(5),
- Hexas[i].x(6), Hexas[i].y(6), Hexas[i].z(6), Hexas[i].x(7), Hexas[i].y(7), Hexas[i].z(7),
- Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag(),Hexas[i].lsTag());
- //Hexas[i].pt(0).ls(),Hexas[i].pt(1).ls(),Hexas[i].pt(2).ls(),Hexas[i].pt(3).ls(),Hexas[i].pt(4).ls(),Hexas[i].pt(5).ls(),Hexas[i].pt(6).ls(),Hexas[i].pt(7).ls());
- }
- int nbIp = ip.size();
- int nbIpIn = 0;
- int nbIpOut = 0;
- for(int i=0;i<nbIp;i++) {
- if(ip[i].isInsideDomain()) nbIpIn++;
- if(ip[i].isOutsideDomain()) nbIpOut++;
- //if(ip[i].isInsideDomain())
- fprintf(f,"SP(%g,%g,%g){%g};\n", ip[i].x(), ip[i].y(), ip[i].z(), ip[i].ls());
- }
-
- int nbIpS = ipS.size();
- for(int i=0;i<nbIpS;i++) {
- //if(ipS[i].x()>=xm && ipS[i].x()<=xM && ipS[i].y()>=ym && ipS[i].y()<=yM && ipS[i].z()>=zm && ipS[i].z()<=zM){
- //if(ipS[i].lsTag()<=0)
- fprintf(f,"SP(%g,%g,%g){%g};\n", ipS[i].x(), ipS[i].y(), ipS[i].z(), ipS[i].ls());//}
- }
- //for(int i=0;i<nbIpS;i++) ipS[i].print();
-
- int nbCP = cp.size();
- for(int i=0;i<nbCP;i++) {
- //fprintf(f,"SP(%g,%g,%g){%g};\n", cp[i].x(), cp[i].y(), cp[i].z(), cp[i].ls());
- }
- fprintf(f,"};\n");
- fclose(f);
-
- printf("%d lines\n", nbLn);
- printf("%d triangles\n", nbTr);
- printf("%d quadrangles\n", nbQ);
- printf("%d tetrahedra\n", nbT);
- printf("%d hexahedra\n", nbH);
-
- printf("%d integration points inside domain\n", nbIpIn);
- printf("%d integration points outside domain\n", nbIpOut);
- printf("%d integration points on boundary\n", nbIpS);
- printf("%d cutting points\n", nbCP);
-
- // Validation by calculating the volume and the surface of a sphere
- double resultVPart=0, resultVTot=0;
- for (int i=0;i<nbIp;i++) {
- resultVTot += ip[i].weight();
- if(ip[i].isInsideDomain()) resultVPart += ip[i].weight();
- }
- double resultSPart=0, resultSTot=0;
- for (int i=0;i<nbIpS;i++) {
- resultSTot += ipS[i].weight();
- if(ipS[i].isOnBorder()) resultSPart += ipS[i].weight();
- }
-
- printf("Integration domain result = %g / %g\n", resultVPart, resultVTot);
- //double volumeS = (4./3.)*(3.141592654*3.45*4.15*4.85)/8.;
- //printf("Exact volume value = %g\n", volumeS);
-
- printf("Integration boundary result = %g / %g\n", resultSPart,resultSTot);
- //double p = 1.6075;
- //double surfaceEllips = 4.0*3.141592654*pow((pow(3.45,p)*pow(4.15,p)+pow(3.45,p)*pow(4.85,p)+pow(4.15,p)*pow(4.85,p))/3.,1./p)/8.;
- //double surfacePart = 5*(4.7836-1.2) + 5*(4.7-0.8) + 5*(3.141592654/2-atan(1.2/4.7)-atan(0.8/4.7836))*4.85;
- //double surfaceTot = 5*5+5*5+5*3.141592654*4.85/2;
- //printf("Exact surface value = %g \n", surfaceEllips);
-
- /*int nMax[4]={1,1,1,1};
- for(int i=0;i<(int)ip.size();i++){
- int nx=digits(ip[i].x());
- int ny=digits(ip[i].y());
- int nz=digits(ip[i].z());
- int nw=digits(ip[i].weight());
- if(nx>nMax[0]) nMax[0]=nx;
- if(ny>nMax[1]) nMax[1]=ny;
- if(nz>nMax[2]) nMax[2]=nz;
- if(nw>nMax[3]) nMax[3]=nw;
- }
-
- printf("nMax=(%d,%d,%d,%d)\n",nMax[0],nMax[1],nMax[2],nMax[3]);
-
- //for(int i=10;i<15;i++) for(int j=10;j<15;j++) printf("%d,%d %-#*.*g\n",i,j,i,j,ip[1].x());
- //for(int i=10;i<15;i++) for(int j=10;j<15;j++) printf("%d,%d %-#*.*g\n",i,j,i,j,ip[2].x());
- for(int i=0;i<(int)ip.size();i++){
- //fprintf(f1," {{%-#12.13g, %-#12.13g, %#.1g}, %-#12.13g},\n",ip[i].x(),ip[i].x(),ip[i].z(),ip[i].x());
- printf(" {{%-#*.*g, %-#*.*g, %-#*.*g}, %+-#*.*g},\n", nMax[0],nMax[0]-nbZeros(ip[i].x()),ip[i].x(), nMax[1],nMax[1]-nbZeros(ip[i].y()),ip[i].y(),
- nMax[2],nMax[2]-nbZeros(ip[i].z()),ip[i].z(), nMax[3],nMax[3]-nbZeros(ip[i].weight()),ip[i].weight());
- }*/
-}
-
-// ----------------------------------------------------------------------
-// ------------------------ READ pMESH ---------------------------------
-// ----------------------------------------------------------------------
-void readMesh(GModel *m, std::vector<DI_Line> &Lines, std::vector<DI_Triangle> &Triangles, std::vector<DI_Quad> &Quads, std::vector<DI_Tetra> &Tetras, std::vector<DI_Hexa> &Hexas)
-{
- // read the mesh and create the tetrahedra and hexahedra or the triangles and quadrangles.
-
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
- for(int i=0;i<(int)(*it)->tetrahedra.size();i++){
- double x[4], y[4], z[4];
- MTetrahedron* t = (*it)->tetrahedra[i];
- for(int j=0;j<4;j++){
- MVertex * v = t->getVertex(j);
- x[j]= v->x();
- y[j]= v->y();
- z[j]= v->z();
- }
- DI_Tetra tt(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3]);
- Tetras.push_back(tt);
- }
- for(int i=0;i<(int)(*it)->hexahedra.size();i++){
- double x[8], y[8], z[8];
- MHexahedron* h = (*it)->hexahedra[i];
- for(int j=0;j<8;j++){
- MVertex * v = h->getVertex(j);
- x[j]= v->x();
- y[j]= v->y();
- z[j]= v->z();
- }
- DI_Hexa hh(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3],x[4],y[4],z[4],x[5],y[5],z[5],x[6],y[6],z[6],x[7],y[7],z[7]); hh.print();
- Hexas.push_back(hh);
- }
- }
-
- /*while (pRegion pr = RIter_next(rit)){
- int numV = R_numVertices(pr);
- double x[numV], y[numV], z[numV];
- pVertex pv;
- pPoint pp;
- for (int i=0; i<numV;++i) {
- pv = R_vertex(pr,i);
- pp = V_point(pv);
- x[i]= P_x(pp);
- y[i]= P_y(pp);
- z[i]= P_z(pp);
- }
- if(numV==4) { //Tetrahedral element
- DI_Tetra tt(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3]);
- Tetras.push_back(tt);
- }
- else if(numV==8){ //Hexahedral element
- DI_Hexa hh(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3],x[4],y[4],z[4],x[5],y[5],z[5],x[6],y[6],z[6],x[7],y[7],z[7]);
- Hexas.push_back(hh);
- }
- else {printf("Volume elements with %d nodes not handled \n",numV); throw;}
- }
- RIter_delete(rit);*/
-
- if(Hexas.size()==0 && Tetras.size()==0) { //add only surfaces if no volume.
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- for(int i=0;i<(int)(*it)->triangles.size();i++){
- double x[3], y[3], z[3];
- MTriangle* t = (*it)->triangles[i];
- for(int j=0;j<3;j++){
- MVertex * v = t->getVertex(j);
- x[j]= v->x();
- y[j]= v->y();
- z[j]= v->z();
- }
- DI_Triangle tt(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2]);
- Triangles.push_back(tt);
- }
- for(int i=0;i<(int)(*it)->quadrangles.size();i++){
- double x[4], y[4], z[4];
- MQuadrangle* q = (*it)->quadrangles[i];
- for(int j=0;j<4;j++){
- MVertex * v = q->getVertex(j);
- x[j]= v->x();
- y[j]= v->y();
- z[j]= v->z();
- }
- DI_Quad qq(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3]);
- Quads.push_back(qq);
- }
- }
- }
-
-
- if(Hexas.size()==0 && Tetras.size()==0 && Triangles.size()==0 && Quads.size()==0) { //add only lines if no surface and no volume.
- for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
- for(int i=0;i<(int)(*it)->lines.size();i++){
- double x[2], y[2], z[2];
- MLine* l = (*it)->lines[i];
- for(int j=0;j<2;j++){
- MVertex * v = l->getVertex(j);
- x[j]= v->x();
- y[j]= v->y();
- z[j]= v->z();
- }
- DI_Line ll(x[0],y[0],z[0],x[1],y[1],z[1]);
- Lines.push_back(ll);
- }
- }
- }
- /*if(Hexas.size()==0 && Tetras.size()==0) { //add only surfaces if no volume.
- Triangles.reserve(M_numTriangles(mesh));
- Quads.reserve(M_numQuads(mesh));
- FIter fit = M_faceIter(mesh);
- while (pFace pf = FIter_next(fit)){
- int numV = F_numVertices(pf);
- double x[numV], y[numV], z[numV];
- pVertex pv;
- pPoint pp;
- for (int i=0; i<numV;++i) {
- pv = F_vertex(pf,i);
- pp = V_point(pv);
- x[i]= P_x(pp);
- y[i]= P_y(pp);
- z[i]= P_z(pp);
- }
- if(numV==3) { // Triangular element
- DI_Triangle tt(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2]);
- Triangles.push_back(tt);
- }
- else if(numV==4) { // Quadrangular element
- DI_Quad qq(x[0],y[0],z[0],x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3]);
- Quads.push_back(qq);
- }
- else {printf("Surface elements with %d nodes not handled \n",numV); throw;}
- }
- FIter_delete(fit);
- }*/
-}
-
-// ----------------------------------------------------------------------
-// -------------MESH INTEGRATION POINTS ---------------------------------
-// ----------------------------------------------------------------------
-
-void getIntegrationPoints(GModel *m, const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip, std::vector<DI_IntegrationPoint> &ipS,
- std::vector<DI_CuttingPoint> &cp, std::vector<DI_Line> &surfLines, std::vector<DI_Triangle> &surfTriangles, std::vector<DI_Quad> &surfQuads,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Hexa> ¬CutHexas)
-{
- clock_t start,end;
- start = clock();
-
- const int POL=1, POTr=1, POQ=1, POT=1, POH=1; //polynomial order for the lines, triangles, the quads, the tetrs and the hexas
-
- std::vector<DI_Line> Lines;
- std::vector<DI_Triangle> Triangles;
- std::vector<DI_Quad> Quads;
- std::vector<DI_Tetra> Tetras;
- std::vector<DI_Hexa> Hexas;
-
-
- readMesh(m, Lines, Triangles, Quads, Tetras, Hexas);
- printf("mesh reading finished : %d Lin, %d Tri, %d Qua, %d Tet, %d Hex\n", Lines.size(), Triangles.size(), Quads.size(), Tetras.size(), Hexas.size());
-
- for(int ln=0;ln<(int)Lines.size();ln++)
- Lines[ln].cut (Ls, ip, cp, POL, surfLines, 0);
-
- for(int tr=0;tr<(int)Triangles.size();tr++)
- Triangles[tr].cut (Ls, ip, ipS, cp, POQ, POTr, POL, surfQuads, surfTriangles, surfLines, 2);
-
- for(int q=0;q<(int)Quads.size();q++) {
- Quads[q].cut (Ls, ip, ipS, cp, POQ, POTr, POL, surfQuads, surfTriangles, surfLines, 2);
- }
-
- for(int t=0;t<(int)Tetras.size();t++){
- Tetras[t].cut (Ls, ip, ipS, cp, POT, POQ, POTr, subTetras, surfQuads, surfTriangles);
- }
-
- for(int h=0;h<(int)Hexas.size();h++){
- Hexas[h].cut (Ls, ip, ipS, cp, POH, POT, POQ, POTr, notCutHexas, subTetras, surfQuads, surfTriangles, surfLines);
- }
-
- end=clock();
- printf("End of computation : %g seconds\n",((double)end - start) / CLOCKS_PER_SEC);
-}
-
-class gLevelsetStrange : public gLevelsetPrimitive
-{
-public:
-
- gLevelsetStrange (int &tag) : gLevelsetPrimitive(tag) {}
- double operator () (const double &x, const double &y, const double &z) const{
- const double X = x-.5;
- const double Y = y-.5;
- const double R0 = .2;
- const double R = sqrt(X*X+Y*Y);
- const double theta = atan2(Y,X);
- const double ls = (R0-R)+(R0/5)*sin(4*theta);
- return ls;
- }
- int type() const {return SPHERE;}
-};
-
-gLevelset *myCircle2b (){
- int tag = 5;
- std::vector<const gLevelset*> li,lj;
- li.push_back(new gLevelsetSphere (.5,.5,0,0.464, tag));
- li.push_back(new gLevelsetStrange(tag));
- lj.push_back(new gLevelsetCut(li));
- lj.push_back(new gLevelsetPlane (1.,0.,0.,-.5,tag));
- gLevelset *ls = new gLevelsetIntersection(lj);
- return ls ;//new gLevelsetReverse(ls);
-}
-
-
-// ---------------------------------------------------------------------
-// ------------------ MAIN FUNCTION ------------------------------------
-// ---------------------------------------------------------------------
-int main(int argc, char **argv)
-{
- int tag=1;
- //double ptc[3]={1,2.5,2.5};
- //double ptc0[3]={3.25,3.25,3.25};
- //double ptc1[3] = {0,1,0};
- //double ptc2[3] = {0,0,1};
- //double ptc3[3] = {2.5,2.5,0.2};
- //double ptc4[3] = {2.5,1.0,2.5};
- //double ptb1[3] = {0,0,0};
- //double ptb2[3] = {1,0,0};
- //double ptb3[3] = {1,1,0};
- //double ptb4[3] = {0,1,0};
- //double ptb5[3] = {0,0,1};
- //double ptb6[3] = {1,0,1};
- //double ptb7[3] = {1,1,1};
- //double ptb8[3] = {0,1,1};
- //double ptr1[3] = {0,-3,0};
- //double dirc[3]={1,0,0};
- //double dir1[3]={1,0,0};
- //double dir2[3]={0,1,0};
- //double dir3[3]={0,0,1};
- //double dir4[3]={-1,0,0};
- //double dir5[3]={0,-1,0};
- //double dir6[3]={0,0,-1};
- const gLevelsetSphere LsS1 (2.5,2.5,0,2.35,tag);
- const gLevelsetSphere LsS2 (2.5,2.5,0,1.75,tag);
- const gLevelsetSphere LsS3 (2.5,5.,0,2.5,tag);
- const gLevelsetPlane LsP1 (0,1,0,-2.9,tag);
- const gLevelsetPlane LsP2 (-1,0,0,1.25,tag);
- //const gLevelsetEllipsoid LsE (ptb1,dir3,2.54,3.6,4.85,tag);
- //const gLevelsetGenCylinder LsCy (ptc,dirc,4.85,tag);
- //const gLevelsetCone LsCo (ptc1,dirc,PI/4,tag);
- std::vector<const gLevelset *> p;
- p.push_back(&LsS1);
- p.push_back(&LsS2);
- //p.push_back(&LsS1);
- //const gLevelsetCut LSC (p);
- p.clear();
- //p.push_back(&LsCy);
- //p.push_back(&LsCo);
- //p.push_back(&LsE);
- //const gLevelsetCut LsCu (p);
- //p.clear();
- //p.push_back(&LsI);
- //p.push_back(&LsCu);
- //const gLevelsetIntersection LsI2 (p);
- //const gLevelset *ls = myCircle2b();
-
- //const gLevelset *ls = new gLevelsetCylinder (ptc,dirc,1.5,1,3,tag);
- //const gLevelsetBox Ls2 (ptc,dir1,dir2,dir3,1.05,0.55,0.8,tag);
- //const gLevelset *ls = new gLevelsetConrod (ptc4, dir2, dir3, 1.4, 0.9, 0.6, 0.8, 0.5, 0.6, 0.3, 1, 0.6, 3,tag);
- //const gLevelsetConrod levset (ptr1,dir2,dir3,2.8,1.8,1.2,1.6,1.,1.2,0.6,2.,1.2,6.,tag);
- //const gLevelsetBox Ls4 (ptb1,ptb2,ptb3,ptb4,ptb5,ptb6,ptb7,ptb8,tag);
- //const gLevelsetCylinder Ls8 (ptc4,dir3,0.8,0.5,1,tag);
-
- //double dirm3[3] = {-dir3[0],-dir3[1],-dir3[2]};
- //double pta2[3] = {ptc3[0],ptc3[1],ptc3[2]+3};
- //std::vector<const gLevelset *> v1;
- //v1.push_back(new gLevelsetGenCylinder (ptc3,dir3,2.2,tag));
- //v1.push_back(new gLevelsetPlane (ptc3,dirm3,tag));
- //v1.push_back(new gLevelsetPlane (pta2,dir3,tag));
- //const gLevelsetIntersection Ls15(v1);
- tag=1;
- //const gLevelset *ls = new gLevelsetPiston(tag);
-
- /*// teeth
- int tag = 5;
- const double p1 [3]={0.125,0.75,0.};
- const double p2 [3]={0.375,0.75,0.};
- const double p3 [3]={0.625,0.75,0.};
- const double p4 [3]={0.875,0.75,0.};
- const double n1 [3]={-1.,0.125/0.75,0.};
- const double n2 [3]={1.,0.125/0.75,0.};
- gLevelset *ls1 = new gLevelsetPlane (p1, n1, tag);
- gLevelset *ls2 = new gLevelsetPlane (p1, n2, tag);
- gLevelset *ls3 = new gLevelsetPlane (p2, n1, tag);
- gLevelset *ls4 = new gLevelsetPlane (p2, n2, tag);
- gLevelset *ls5 = new gLevelsetPlane (p3, n1, tag);
- gLevelset *ls6 = new gLevelsetPlane (p3, n2, tag);
- gLevelset *ls7 = new gLevelsetPlane (p4, n1, tag);
- gLevelset *ls8 = new gLevelsetPlane (p4, n2, tag);
- std::vector<const gLevelset*> l1,l2,l3,l4,li;
- l1.push_back(ls1);
- l1.push_back(ls2);
- l2.push_back(ls3);
- l2.push_back(ls4);
- l3.push_back(ls5);
- l3.push_back(ls6);
- l4.push_back(ls7);
- l4.push_back(ls8);
- li.push_back(new gLevelsetIntersection(l1));
- li.push_back(new gLevelsetIntersection(l2));
- li.push_back(new gLevelsetIntersection(l3));
- li.push_back(new gLevelsetIntersection(l4));
- const gLevelset *ls = new gLevelsetUnion(li);*/
-
- /*const double p1 [3]={1.,2.5,2.5};
- const double n1 [3]={1.,0.,0.};
- const double n2 [3]={0.,1.,0.};
- const double n3 [3]={0.,0.,1.};
- gLevelset *ls1 = new gLevelsetPlane (p1, n1, tag);
- gLevelset *ls2 = new gLevelsetPlane (p1, n2, tag);
- gLevelset *ls3 = new gLevelsetPlane (p1, n3, tag);
- std::vector<const gLevelset*> l1;
- l1.push_back(ls1);
- l1.push_back(ls2);
- l1.push_back(ls3);
- const gLevelset *ls = new gLevelsetUnion(l1);*/
-
- const double pt [3]={0.,0.,0.};
- const double dir[3]={0,0,1};
- tag=1;
- std::vector<const gLevelset*> li;
- li.push_back(new gLevelsetPlane(0,-1,0,0.52, tag)); // -/+
- li.push_back(new gLevelsetPlane(-1,0,0,0.52, tag));
- //li.push_back(new gLevelsetGenCylinder(pt,dir,.75, tag)); // -|+
- gLevelset *ls = new gLevelsetIntersection(li);
- //const gLevelset *ls = new gLevelsetReverse (ls1);
- //gLevelset *ls = new gLevelsetGenCylinder(pt, dir, 0.75, tag);
-
- GmshInitialize(argc, argv);
- GmshSetOption("General", "Terminal", 1.);
- GModel *m = new GModel;
- m->readMSH(argv[1]); printf("mesh read\n");
-
- /*printf(" carreTri\n");
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- printf("ptf %d (%g,%g)\n",(*it)->mesh_vertices[i]->getNum(),(*it)->mesh_vertices[i]->x(),(*it)->mesh_vertices[i]->y());
- for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- printf("pte %d (%g,%g)\n",(*it)->mesh_vertices[i]->getNum(),(*it)->mesh_vertices[i]->x(),(*it)->mesh_vertices[i]->y()); printf("\n");*/
-
- GModel *cm = m->buildCutGModel(ls); printf("mesh cut\n");
-
- /*printf(" ccutmesh\n");
- for(GModel::fiter it = cm->firstFace(); it != cm->lastFace(); ++it){
- printf("face %d\n",(*it)->tag());
- for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
- printf("tri %d (%d, %d, %d)\n",(*it)->triangles[i]->getNum(),(*it)->triangles[i]->getVertex(0)->getNum(),(*it)->triangles[i]->getVertex(1)->getNum(),(*it)->triangles[i]->getVertex(2)->getNum());
- //for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
- //printf("poly %d (%d, %d, %d)\n",(*it)->polygons[i]->getNum(),(*it)->polygons[i]->getVertex(0)->getNum(),(*it)->polygons[i]->getVertex(1)->getNum(),(*it)->polygons[i]->getVertex(2)->getNum());
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- printf("ptf %d (%g,%g)\n",(*it)->mesh_vertices[i]->getNum(),(*it)->mesh_vertices[i]->x(),(*it)->mesh_vertices[i]->y());
- }
- for(GModel::eiter it = cm->firstEdge(); it != cm->lastEdge(); ++it){
- printf("edge %d\n",(*it)->tag());
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- printf("pte %d (%g,%g)\n",(*it)->mesh_vertices[i]->getNum(),(*it)->mesh_vertices[i]->x(),(*it)->mesh_vertices[i]->y());
- }*/
-
- cm->writeMSH("cutMesh.msh",2); printf("mesh written\n");
-
- /*GFace *gf = m->getFaceByTag(6); std::cout<<"face "<<gf<< " "<<gf->mesh_vertices.size()<<std::endl;
- GEdge *ge1 = m->getEdgeByTag(1); std::cout<<"edge1 "<<ge1<< " "<<ge1->mesh_vertices.size()<<std::endl;
- GEdge *ge2 = m->getEdgeByTag(2); std::cout<<"edge2 "<<ge2<< " "<<ge2->mesh_vertices.size()<<std::endl;
- GEdge *ge3 = m->getEdgeByTag(3); std::cout<<"edge3 "<<ge3<< " "<<ge3->mesh_vertices.size()<<std::endl;
- GEdge *ge4 = m->getEdgeByTag(4); std::cout<<"edge4 "<<ge4<< " "<<ge4->mesh_vertices.size()<<std::endl;
- for(int i = 0; i < gf->mesh_vertices.size(); i++)
- if(gf->mesh_vertices[i]->x() == 2 && gf->mesh_vertices[i]->y() == 0)
- std::cout << "face MVertex " << gf->mesh_vertices[i]->getNum() << " " << gf->mesh_vertices[i]->onWhat()->tag() << " " << gf->mesh_vertices[i] << " " << gf->mesh_vertices[i]->onWhat() << std::endl;
- for(int i = 0; i < ge3->mesh_vertices.size(); i++)
- if(ge3->mesh_vertices[i]->x() == 2 && ge3->mesh_vertices[i]->y() == 0)
- std::cout << "edge3 MVertex " << ge3->mesh_vertices[i]->getNum() << " " << ge3->mesh_vertices[i]->onWhat()->tag() << " " << ge3->mesh_vertices[i] << " " << ge3->mesh_vertices[i]->onWhat() << std::endl;
- for(int i = 0; i < ge4->mesh_vertices.size(); i++)
- if(ge4->mesh_vertices[i]->x() == 2 && ge4->mesh_vertices[i]->y() == 0)
- std::cout << "edge4 MVertex " << ge4->mesh_vertices[i]->getNum() << " " << ge4->mesh_vertices[i]->onWhat()->tag() << " " << ge4->mesh_vertices[i] << " " << ge4->mesh_vertices[i]->onWhat() << std::endl;*/
-
- //pGModel model = 0;
- //pMesh mesh = MS_newMesh(model);
- //M_load(mesh, "cubeHex.msh");
-
- std::vector<DI_IntegrationPoint> ip; // contains the volume Integration points
- std::vector<DI_IntegrationPoint> ipS; // contains the surface Integration points
- std::vector<DI_CuttingPoint> cp; // contains the cutting points
- std::vector<DI_Line> surfLines; // contains the line contour
- std::vector<DI_Triangle> surfTriangles; // contains the triangles forming the crack surface
- std::vector<DI_Quad> surfQuads; // contains the quads forming the crack surface
- std::vector<DI_Tetra> Tetras; // contains the tetrahedra formed by the cutting of the hexahedra/tetrahedra and the tetrahedra not cut by the level set
- std::vector<DI_Hexa> Hexas; // contains the hexahedra not cut by the level set
-
- std::vector<const gLevelset *> primitives;
- ls->getPrimitives(primitives);
- //printf("prim = \n"); for(int i=0;i<(int)primitives.size();i++) {primitives[i]->print(); printf("LS(0,0)=%g LS(1,1)=%g\n",(*primitives[i])(0,0,0),(*primitives[i])(1,1,0));}
- std::vector<const gLevelset *> RPN;
- ls->getRPN(RPN);
- //printf("RPN = \n"); for(int i=0;i<(int)RPN.size();i++) {RPN[i]->print(); printf("LS(0,0)=%g LS(1,1)=%g\n",(*RPN[i])(0,0,0),(*RPN[i])(1,1,0));}
-
-// getIntegrationPoints(m, *ls, ip, ipS, cp, surfLines, surfTriangles, surfQuads, Tetras, Hexas);
-// writePosFile (ip, ipS, cp, surfLines, surfTriangles, surfQuads, Tetras, Hexas);
-
-}