Galerkin FEM for elliptic PDEs
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

sphere.geo 4.2KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. radius = 3.25; // Radius of the damn thing
  9. lc = radius/5; // 0.25; // Target element size
  10. // Total Solution Space
  11. Box = 3*radius; // The down side of potential
  12. X0 = -Box;
  13. X1 = Box;
  14. Y0 = -Box;
  15. Y1 = Box;
  16. Z0 = -Box;
  17. Z1 = Box;
  18. cellSize=radius/10; ///10;
  19. dd = 0 ; // 1e-5; //cellSize; // .01;
  20. pio2=Pi/2;
  21. /////////////////////////////////////
  22. // Large Bounding box
  23. pp = newp;
  24. Point(pp) = {X0, Y0, Z0, lc};
  25. Point(pp+1) = {X1, Y0, Z0, lc};
  26. Point(pp+2) = {X1, Y1, Z0, lc};
  27. Point(pp+3) = {X0, Y1, Z0, lc};
  28. lv = newl;
  29. Line(lv) = {pp,pp+1};
  30. Line(lv+1) = {pp+1,pp+2};
  31. Line(lv+2) = {pp+2,pp+3};
  32. Line(lv+3) = {pp+3,pp};
  33. Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
  34. // Hard coded doom
  35. Plane Surface(125) = {lv+4};
  36. //v = newv;
  37. v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
  38. // Calculate offset effect
  39. theta = Asin(dd/radius);
  40. rr = radius * Cos(theta);
  41. ///////////////////////////////////
  42. // Positive half sphere
  43. // create inner 1/8 shell
  44. p0 = newp;
  45. Point(p0) = { 0, 0, 0, cellSize}; // origin
  46. Point(p0+1) = { -rr, 0, dd, cellSize};
  47. Point(p0+2) = { 0, rr, dd, cellSize};
  48. Point(p0+3) = { 0, 0, radius, cellSize};
  49. Point(p0+4) = { 0, 0, dd, cellSize}; // origin
  50. c0 = newc;
  51. Circle(c0 ) = {p0+1, p0+4, p0+2}; // Tricky, This one needs to be offset!
  52. Circle(c0+1) = {p0+3, p0, p0+1};
  53. Circle(c0+2) = {p0+3, p0, p0+2};
  54. Line Loop(10) = {c0, -(c0+2), c0+1} ;
  55. Ruled Surface (60) = {10};
  56. ////////////////////////////////////////////////////////////
  57. // Negative half sphere
  58. p = newp;
  59. Point( p) = { 0, 0, 0, cellSize};
  60. Point(p+1) = { -rr, 0, -dd, cellSize};
  61. Point(p+2) = { 0, rr, -dd, cellSize};
  62. Point(p+3) = { 0, 0, -radius, cellSize};
  63. Point(p+4) = { 0, 0, -dd, cellSize};
  64. cc = newc;
  65. Circle(cc ) = {p+1, p+4, p+2};
  66. Circle(cc+1) = {p+3, p, p+1};
  67. Circle(cc+2) = {p+3, p, p+2};
  68. Circle(cc+3) = {p+3, p, p+2};
  69. Circle(cc+4) = {p+3, p, p+2};
  70. Circle(cc+5) = {p+3, p, p+2};
  71. ccl = newl;
  72. Line(ccl) = { p0+3, p+3 };
  73. Line Loop(11) = {cc, -(cc+2), cc+1} ;
  74. Ruled Surface (61) = {11};
  75. // create remaining 7/8 inner shells
  76. t1[] = Rotate {{0,0,1},{0,0,0},pio2 } {Duplicata{Surface{60};}};
  77. t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
  78. t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
  79. //
  80. t4[] = Rotate {{0,0,1},{0,0,0},pio2 } {Duplicata{Surface{61};}};
  81. t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
  82. t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
  83. /* This is GOOD */
  84. Surface{60} In Volume{v[1]};
  85. Surface{t1[0]} In Volume{v[1]};
  86. Surface{t2[0]} In Volume{v[1]};
  87. Surface{t3[0]} In Volume{v[1]};
  88. Surface{61} In Volume{v[1]};
  89. Surface{t4[0]} In Volume{v[1]};
  90. Surface{t5[0]} In Volume{v[1]};
  91. Surface{t6[0]} In Volume{v[1]};
  92. ///////////////////////////////////////////////
  93. // Attractor Field
  94. // Field[1] = Attractor;
  95. // Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
  96. //
  97. // Field[2] = MathEval;
  98. // Field[2].F = Sprintf("(2.25 - F1)^1.01 + %g", radius/10 ); // WORKS
  99. //
  100. // Field[3] = MathEval;
  101. // Field[3].F = Sprintf("(12.25 - F1)^1.01 + %g", radius/5 ); // WORKS
  102. //
  103. // Field[4] = MathEval;
  104. // Field[4].F = Sprintf("(22.25 - F1)^1.01 + %g", radius ); // WORKS
  105. //
  106. // Field[5] = MathEval;
  107. // Field[5].F = Sprintf("(42.25 - F1)^1.01 + %g", radius*5 ); // WORKS
  108. //
  109. // //Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, 2*cellSize );
  110. // //Background Field = 2;
  111. //
  112. // // Finally, let's use the minimum of all the fields as the background mesh field
  113. // Field[7] = Min;
  114. // Field[7].FieldsList = {2, 3, 4, 5};
  115. // Background Field = 7;
  116. // Don't extend the elements sizes from the boundary inside the domain
  117. //Mesh.CharacteristicLengthExtendFromBoundary = 0;
  118. Physical Volume(1) = {v[1]};
  119. // To create the mesh run
  120. // gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
  121. //gmsh -3 -format msh1 -o outfile.msh sphere.geo