[Gmsh] Suspect self-intersecting curvilinear geometry
Aleksejs Fomins
aleksejs.fomins at lspr.ch
Wed Jun 17 16:09:12 CEST 2015
Dear GMSH,
I have been peacefully using GMSH-generated curvilinear geometries (up to order 5) for some time now using my hand-written GMSH reader. It seemed to work fine on most curvilinear geometries I have generated so far.
However, recently I managed to generate a 2nd order curvilinear geometry, for which my self-test fails. The self-test claims that the jacobian determinant of one of the elements changes sign over the element, which is a clear indication that the element is self-intersecting.
>From what I have understood from the documentation, this case should be impossible.
I will attach the generated .msh and .geo files. It would be great if somebody could confirm if the mesh indeed possesses self-intersecting elements, and gave me suggestions on what could have gone wrong.
Kind regards,
Aleksejs Fomins
-------------- next part --------------
/** \file
* \brief Gmsh script that generates three nested spheres with outermost radius=1 unit
*
* Copyright by Patrick Leidenberger and Benedikt Oswald, 2006-2009.
* All rights reserved.
*
* Objective: Gmsh script of a sphere in a sphere in a sphere.
*
* \author Patrick Leidenberger, Benedikt Oswald
* \date 2006 dec 15, created, Patrick Leidenberger.
* \date 2006 dec 15, modified, Patrick Leidenberger,
* \date 2007 apr 26, modified, Benedikt Oswald, added physical boundary and volume tags
* \date 2009 aug 04, modified, Benedikt Oswald, changes to boundary and volume tags
*
* \warning None.
* \attention
* \bug
* \todo
../../hadesgeo/babel/gmsh2dgf/gmsh2dgf \
--input-file=triple-embedded-sphere-curvilinear-revision-1000.msh \
--output-file=triple-embedded-sphere-curvilinear-revision-1000.dgf \
--tag-original-range-start=100 --tag-original-range-stop=150 --tag-target-value=101 \
--tag-original-range-start=500 --tag-original-range-stop=900 --tag-target-value=501 \
--tag-original-range-start=901 --tag-original-range-stop=920 --tag-target-value=901 \
--tag-original-range-start=950 --tag-original-range-stop=990 --tag-target-value=950
*/
// define physical tags
innermost_surface = 301;
inbetween_surface = 201;
outermost_surface = 101;
innermost_volume = 501;
inbetween_volume = 901;
outermost_volume = 950;
scaling = 8.0; // Scaling factor for all.
cl = scaling * 1.00; // Characteristic length for all.
cl0 = scaling * 1.00; // Characteristic length sphere 0 = inner most sphere
cl1 = scaling * 1.00; // Characteristic length sphere 1 = sphere in between inner most and outer most sphere
cl2 = scaling * 1.00; // Characteristic length sphere 2 = outer most sphere
// Center of sphere.
centerX = 0.0;
centerY = 0.0;
centerZ = 0.0;
// Radius of the sphere
radius0 = 0.050; // radius of inner most sphere
radius1 = 0.165; // radius of sphere in between inner most and outer most sphere
radius2 = 0.200; // radius of outer most sphere
// Create points on sphere 0 surface = innermost sphere
ipt1 = newp; Point(ipt1) = {centerX, centerY, centerZ, cl0};
ipt2 = newp; Point(ipt2) = {centerX - radius0, centerY, centerZ, cl0};
ipt3 = newp; Point(ipt3) = {centerX + radius0, centerY, centerZ, cl0};
ipt4 = newp; Point(ipt4) = {centerX, centerY - radius0, centerZ, cl0};
ipt5 = newp; Point(ipt5) = {centerX, centerY + radius0, centerZ, cl0};
ipt6 = newp; Point(ipt6) = {centerX, centerY, centerZ - radius0, cl0};
ipt7 = newp; Point(ipt7) = {centerX, centerY, centerZ + radius0, cl0};
// Create points on sphere 1 surface = spherical shell in between
ipt8 = newp; Point(ipt8) = {centerX, centerY, centerZ, cl1};
ipt9 = newp; Point(ipt9) = {centerX - radius1, centerY, centerZ, cl1};
ipt10 = newp; Point(ipt10) = {centerX + radius1, centerY, centerZ, cl1};
ipt11 = newp; Point(ipt11) = {centerX, centerY - radius1, centerZ, cl1};
ipt12 = newp; Point(ipt12) = {centerX, centerY + radius1, centerZ, cl1};
ipt13 = newp; Point(ipt13) = {centerX, centerY, centerZ - radius1, cl1};
ipt14 = newp; Point(ipt14) = {centerX, centerY, centerZ + radius1, cl1};
// Create points on sphere 2 surface = outermost spherical shell
ipt15 = newp; Point(ipt15) = {centerX, centerY, centerZ, cl2};
ipt16 = newp; Point(ipt16) = {centerX - radius2, centerY, centerZ, cl2};
ipt17 = newp; Point(ipt17) = {centerX + radius2, centerY, centerZ, cl2};
ipt18 = newp; Point(ipt18) = {centerX, centerY - radius2, centerZ, cl2};
ipt19 = newp; Point(ipt19) = {centerX, centerY + radius2, centerZ, cl2};
ipt20 = newp; Point(ipt20) = {centerX, centerY, centerZ - radius2, cl2};
ipt21 = newp; Point(ipt21) = {centerX, centerY, centerZ + radius2, cl2};
// Create circle sections connecting two points on sphere 0 surface.
icl1 = newreg; Circle(icl1) = {ipt2,ipt1,ipt4};
icl2 = newreg; Circle(icl2) = {ipt2,ipt1,ipt5};
icl3 = newreg; Circle(icl3) = {ipt2,ipt1,ipt6};
icl4 = newreg; Circle(icl4) = {ipt2,ipt1,ipt7};
icl5 = newreg; Circle(icl5) = {ipt3,ipt1,ipt4};
icl6 = newreg; Circle(icl6) = {ipt3,ipt1,ipt5};
icl7 = newreg; Circle(icl7) = {ipt3,ipt1,ipt6};
icl8 = newreg; Circle(icl8) = {ipt3,ipt1,ipt7};
icl9 = newreg; Circle(icl9) = {ipt4,ipt1,ipt6};
icl10 = newreg; Circle(icl10) = {ipt4,ipt1,ipt7};
icl11 = newreg; Circle(icl11) = {ipt5,ipt1,ipt6};
icl12 = newreg; Circle(icl12) = {ipt5,ipt1,ipt7};
// Create circle sections connecting two points on sphere 1 surface.
icl13 = newreg; Circle(icl13) = {ipt9,ipt8,ipt11};
icl14 = newreg; Circle(icl14) = {ipt9,ipt8,ipt12};
icl15 = newreg; Circle(icl15) = {ipt9,ipt8,ipt13};
icl16 = newreg; Circle(icl16) = {ipt9,ipt8,ipt14};
icl17 = newreg; Circle(icl17) = {ipt10,ipt8,ipt11};
icl18 = newreg; Circle(icl18) = {ipt10,ipt8,ipt12};
icl19 = newreg; Circle(icl19) = {ipt10,ipt8,ipt13};
icl20 = newreg; Circle(icl20) = {ipt10,ipt8,ipt14};
icl21 = newreg; Circle(icl21) = {ipt11,ipt8,ipt13};
icl22 = newreg; Circle(icl22) = {ipt11,ipt8,ipt14};
icl23 = newreg; Circle(icl23) = {ipt12,ipt8,ipt13};
icl24 = newreg; Circle(icl24) = {ipt12,ipt8,ipt14};
// Create circle sections connecting two points on sphere 2 surface.
icl25 = newreg; Circle(icl25) = {ipt16,ipt15,ipt18};
icl26 = newreg; Circle(icl26) = {ipt16,ipt15,ipt19};
icl27 = newreg; Circle(icl27) = {ipt16,ipt15,ipt20};
icl28 = newreg; Circle(icl28) = {ipt16,ipt15,ipt21};
icl29 = newreg; Circle(icl29) = {ipt17,ipt15,ipt18};
icl30 = newreg; Circle(icl30) = {ipt17,ipt15,ipt19};
icl31 = newreg; Circle(icl31) = {ipt17,ipt15,ipt20};
icl32 = newreg; Circle(icl32) = {ipt17,ipt15,ipt21};
icl33 = newreg; Circle(icl33) = {ipt18,ipt15,ipt20};
icl34 = newreg; Circle(icl34) = {ipt18,ipt15,ipt21};
icl35 = newreg; Circle(icl35) = {ipt19,ipt15,ipt20};
icl36 = newreg; Circle(icl36) = {ipt19,ipt15,ipt21};
// Make the surface mesh for 3 closed circle sections; sphere0.
ill1 = newreg; Line Loop(ill1) = {icl1,-icl3,icl9} ; irs1= newreg; Ruled Surface(irs1) = {ill1};
ill2 = newreg; Line Loop(ill2) = {icl1,icl10,-icl4} ; irs2= newreg; Ruled Surface(irs2) = {ill2};
ill3 = newreg; Line Loop(ill3) = {icl2,-icl3,icl11} ; irs3= newreg; Ruled Surface(irs3) = {ill3};
ill4 = newreg; Line Loop(ill4) = {icl2,icl12,-icl4} ; irs4= newreg; Ruled Surface(irs4) = {ill4};
ill5 = newreg; Line Loop(ill5) = {icl5,icl9,-icl7} ; irs5= newreg; Ruled Surface(irs5) = {ill5};
ill6 = newreg; Line Loop(ill6) = {icl5,icl10,-icl8} ; irs6= newreg; Ruled Surface(irs6) = {ill6};
ill7 = newreg; Line Loop(ill7) = {icl6,icl11,-icl7} ; irs7= newreg; Ruled Surface(irs7) = {ill7};
ill8 = newreg; Line Loop(ill8) = {icl6,icl12,-icl8} ; irs8= newreg; Ruled Surface(irs8) = {ill8};
// Make the surface mesh for 3 closed circle sections; sphere1.
ill9 = newreg; Line Loop(ill9) = {icl13,-icl15,icl21} ; irs9= newreg; Ruled Surface(irs9) = {ill9};
ill10 = newreg; Line Loop(ill10) = {icl13,icl22,-icl16} ; irs10= newreg; Ruled Surface(irs10) = {ill10};
ill11 = newreg; Line Loop(ill11) = {icl14,-icl15,icl23} ; irs11= newreg; Ruled Surface(irs11) = {ill11};
ill12 = newreg; Line Loop(ill12) = {icl14,icl24,-icl16} ; irs12= newreg; Ruled Surface(irs12) = {ill12};
ill13 = newreg; Line Loop(ill13) = {icl17,icl21,-icl19} ; irs13= newreg; Ruled Surface(irs13) = {ill13};
ill14 = newreg; Line Loop(ill14) = {icl17,icl22,-icl20} ; irs14= newreg; Ruled Surface(irs14) = {ill14};
ill15 = newreg; Line Loop(ill15) = {icl18,icl23,-icl19} ; irs15= newreg; Ruled Surface(irs15) = {ill15};
ill16 = newreg; Line Loop(ill16) = {icl18,icl24,-icl20} ; irs16= newreg; Ruled Surface(irs16) = {ill16};
// Make the surface mesh for 3 closed circle sections; sphere2.
ill17 = newreg; Line Loop(ill17) = {icl25,-icl27,icl33} ; irs17= newreg; Ruled Surface(irs17) = {ill17};
ill18 = newreg; Line Loop(ill18) = {icl25,icl34,-icl28} ; irs18= newreg; Ruled Surface(irs18) = {ill18};
ill19 = newreg; Line Loop(ill19) = {icl26,-icl27,icl35} ; irs19= newreg; Ruled Surface(irs19) = {ill19};
ill20 = newreg; Line Loop(ill20) = {icl26,icl36,-icl28} ; irs20= newreg; Ruled Surface(irs20) = {ill20};
ill21 = newreg; Line Loop(ill21) = {icl29,icl33,-icl31} ; irs21= newreg; Ruled Surface(irs21) = {ill21};
ill22 = newreg; Line Loop(ill22) = {icl29,icl34,-icl32} ; irs22= newreg; Ruled Surface(irs22) = {ill22};
ill23 = newreg; Line Loop(ill23) = {icl30,icl35,-icl31} ; irs23= newreg; Ruled Surface(irs23) = {ill23};
ill24 = newreg; Line Loop(ill24) = {icl30,icl36,-icl32} ; irs24= newreg; Ruled Surface(irs24) = {ill24};
ilsl1 = newreg;
Surface Loop(ilsl1) = {ill1+1,ill2+1,ill3+1,ill4+1,ill5+1,ill6+1,ill7+1,ill8+1};
ilsl2 = newreg;
Surface Loop(ilsl2) = {ill1+1,ill2+1,ill3+1,ill4+1,ill5+1,ill6+1,ill7+1,ill8+1,ill9+1,ill10+1,ill11+1,ill12+1,ill13+1,ill14+1,ill15+1,ill16+1};
ilsl3 = newreg;
Surface Loop(ilsl3) = {ill9+1,ill10+1,ill11+1,ill12+1,ill13+1,ill14+1,ill15+1,ill16+1,ill17+1,ill18+1,ill19+1,ill20+1,ill21+1,ill22+1,ill23+1,ill24+1};
// define physical surface tags
//Physical Surface (innermost_surface) = {irs1,irs2,irs3,irs4,irs5,irs6,irs7,irs8};
//Physical Surface (inbetween_surface) = {irs9,irs10,irs11,irs12,irs13,irs14,irs15,irs16};
Physical Surface (outermost_surface) = {irs17,irs18,irs19,irs20,irs21,irs22,irs23,irs24};
// define volume of sphere.
ivl1 = newv; Volume(ivl1) = {ilsl1};
ivl2 = newv; Volume(ivl2) = {ilsl2};
ivl3 = newv; Volume(ivl3) = {ilsl3};
// define physical volumes
Physical Volume (innermost_volume) = ivl1;
Physical Volume (inbetween_volume) = ivl2;
Physical Volume (outermost_volume) = ivl3;
Coherence;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: triple-embedded-sphere-curvilinear-revision-1010-order2.msh
Type: model/mesh
Size: 49999 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20150617/e0790e76/attachment.msh>