[Gmsh] projecting points to geometry
Rene Schneider
rene.schneider at mathematik.tu-chemnitz.de
Thu Mar 21 14:23:13 CET 2013
Dear list
I'm developing a FEM solver, and recently integrated a routine to read
gmsh .msh files. The FEM solver uses hierarchical meshes to enable fast
solvers (multigrid preconditioning). For this reason I prefer to start
with a coarse mesh (generated by gmsh) and to let my code refine this mesh.
In order to deal with curved geometry in the refinement, it is common
practise to project the new surface points onto their (curved) geometry.
Since I didn't find a simple means to do so with gmsh directly, I
modified one of the api-demos to this end. See attached files (which I
keep in the utils/api_demos directory).
This routine reads a file with a list of points and associated geometry
tags. Then it uses the
edge->closestPoint
and
face->closestPoint
methods, to project these points onto their respective geometry. The
result is then written to file again, for the FE solver to read.
So far this seems to be working fine. (See below for a performance example.)
A few comments/questions:
1. I'd be glad if some gmsh developers could comment on the code.
(Is there a better way to achieve this?)
2. My attempts in the CMakeList.txt to create a statically linked
version all failed. Any suggestions on this?
3. I provide these files under same license as gmsh in the hope that
this can be integrated into future releases of gmsh. I suggest the
utils/api_demos folder.
4. Is it easily possible to influence the accuracy of the projection? I
had a look on some of the code for closestPoint and noticed that this is
done in an iterative way, with fairly loose stoping criterion. In the
future I might be interested in differentiating the node positions with
respect to geometry parameters. These derivatives could easily be
approximated by finite differences. However, the projection accuracy is
critical then.
5. Are there some error-flags in case the projection goes wrong? For one
geometry
http://www.opencascade.org/ex/att/15_cylinder_head.brep.gz
from
http://www.opencascade.org/showroom/shapegallery/gal4/
I had problems with a some of the sub-geometries, which were projected
to the origin. Even though I got no error messages or similar.
Regards
Rene
P.S.: As a sample of the performance:
Poisson equation in 3D domain with curved boundary from
http://www.opencascade.org/ex/att/31_misc2.brep.gz
We use gmsh to create a coarse mesh of 7,334 nodes, P1 elements. The
finest mesh so far has 21,909,034 nodes. Solution on this finest mesh is
computed in 15m34.028s including file write of solution for paraview. :)
Visualisation (paraview):
http://www-user.tu-chemnitz.de/~rens/software/feins/examples/full/31_misc2_lvl4_21909034_nodes_s.png
http://www-user.tu-chemnitz.de/~rens/software/feins/examples/full/31_misc2_lvl4_21909034_nodes_swe_large.png
--
----------------------------------------------------
Dr. Rene Schneider
TU Chemnitz, Fakultaet fuer Mathematik,
09107 Chemnitz, Germany
Besucheradresse / Visitor address:
Reichenhainer Str. 41 / Raum 625
09126 Chemnitz, Germany
Tel.: +49-371-531-33953
Fax: +49-371-531-8-33953
rene.schneider at mathematik.tu-chemnitz.de
http://www.tu-chemnitz.de/~rens
----------------------------------------------------
-------------- next part --------------
#include <stdio.h>
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
int main(int argc, char **argv)
{
GmshInitialize(argc, argv);
if (argc<5)
{
fprintf(stderr,"\n%s: needs at least 4 arguments, "
"usage\n %s name.geo name.msh querry.qrx output.pnt\n", argv[0], argv[0]);
exit(1);
}
GModel *m = new GModel();
m->readGEO(argv[1]);
m->readMSH(argv[2]);
/* open querry file for reading */
FILE *querry = fopen(argv[3], "r");
if (querry == NULL)
{
fprintf(stderr, "%s: couldn't open querry file: %s\n",
argv[0],argv[3]);
exit(1);
}
FILE *output = fopen(argv[4], "w");
if (output == NULL)
{
fprintf(stderr, "%s: couldn't open output file: %s\n",
argv[0],argv[4]);
exit(1);
}
/* expected querry format:
k-number-line-points
index1 line-tag current-x current-y current-z
...
indexk line-tag current-x current-y current-z
m-number-face-points
index1 face-tag current-x current-y current-z
...
indexm face-tag current-x current-y current-z
[EOF]
output format is the same, but without face-tag
*/
int i, numEdgePoints, numFacePoints, nnread;
nnread = fscanf(querry,"%d", &numEdgePoints);
if (nnread!=1)
{
fprintf(stderr, "%s: error reading number of edges in file: %s\n",
argv[0], argv[3]);
exit(1);
}
fprintf(output,"%d\n", numEdgePoints);
for (i=0; i<numEdgePoints; i++)
{
int id, lineTag, nread;
double qx,qy,qz;
nread = fscanf(querry,"%d %d %le %le %le ", &id, &lineTag, &qx, &qy, &qz);
if (nread!=5)
{
fprintf(stderr, "%s: error reading %d-th edge entry in file: %s\n",
argv[0], i+1, argv[3]);
exit(1);
}
GEdge *eg = m->getEdgeByTag(lineTag);
SPoint3 qpoint(qx,qy,qz);
double param;
GPoint closest = eg->closestPoint(qpoint, param);
fprintf(output,"%d %25.16e %25.16e %25.16e\n",
id, closest.x(), closest.y(), closest.z());
}
nnread = fscanf(querry,"%d", &numFacePoints);
if (nnread!=1)
{
fprintf(stderr, "%s: error reading number of face in file: %s\n",
argv[0], argv[3]);
exit(1);
}
fprintf(output,"%d\n", numFacePoints);
for (i=0; i<numFacePoints; i++)
{
int id, faceTag, nread;
double qx,qy,qz;
nread = fscanf(querry,"%d %d %le %le %le ", &id, &faceTag, &qx, &qy, &qz);
if (nread!=5)
{
fprintf(stderr, "%s: error reading %d-th face entry in file: %s\n",
argv[0], i+1, argv[3]);
exit(1);
}
GFace *fc = m->getFaceByTag(faceTag);
SPoint3 qpoint(qx,qy,qz);
double initialguess[2] = {0.0, 0.0};
GPoint closest = fc->closestPoint(qpoint, initialguess);
fprintf(output,"%d %25.16e %25.16e %25.16e\n",
id, closest.x(), closest.y(), closest.z());
}
delete m;
GmshFinalize();
exit(0);
}
-------------- next part --------------
cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
# if CMAKE_BUILD_TYPE is specified use it; otherwise set the default
# build type to "RelWithDebInfo" ("-O2 -g" with gcc) prior to calling
# project()
if(DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE} CACHE STRING "Choose build type")
else(DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE RelWithDebInfo CACHE STRING "Choose build type")
endif(DEFINED CMAKE_BUILD_TYPE)
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}
${CMAKE_CURRENT_BINARY_DIR}/gmsh/Common)
if(APPLE)
set(glut "-framework GLUT")
else(APPLE)
set(glut "glut")
endif(APPLE)
add_executable(mainAntTweakBar mainAntTweakBar.cpp)
target_link_libraries(mainAntTweakBar shared AntTweakBar ${glut})
add_executable(mainCartesian mainCartesian.cpp)
target_link_libraries(mainCartesian shared)
add_executable(mainElasticity mainElasticity.cpp)
target_link_libraries(mainElasticity shared)
add_executable(mainGlut mainGlut.cpp)
target_link_libraries(mainGlut lib ${GMSH_EXTERNAL_LIBRARIES} ${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)
add_executable(mainPost mainPost.cpp)
target_link_libraries(mainPost shared)
add_executable(mainSimple mainSimple.cpp)
target_link_libraries(mainSimple shared)
add_executable(mainFEINStest mainFEINStest.cpp)
target_link_libraries(mainFEINStest shared)
add_executable(gmsh_closest_points_shared gmsh_closest_points.cpp)
target_link_libraries(gmsh_closest_points_shared shared)
add_executable(gmsh_closest_points_static gmsh_closest_points.cpp)
target_link_libraries(gmsh_closest_points_static lib ${GMSH_EXTERNAL_LIBRARIES} ${LINK_LIBRARIES})
message("gmsh_src: " ${GMSH_SRC})
add_executable(gmsh_closest_points_src gmsh_closest_points.cpp ${GMSH_SRC})
target_link_libraries(gmsh_closest_points_src ${LINK_LIBRARIES})