Galerkin FEM for elliptic PDEs
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

toroid.geo 3.0KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105
  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. /**
  9. * @file
  10. * @date 02/04/2016 02:58:54 PM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email tirons@egi.utah.edu
  14. * @copyright Copyright (c) 2016, University of Utah
  15. * @copyright Copyright (c) 2016, Lemma Software, LLC
  16. */
  17. radius = 3.25; // Radius of the damn thing
  18. blc = radius/2; // 0.25; // Target element size
  19. Box = 3*radius; // The down side of potential
  20. lc = radius/2; // toroid characteristic length
  21. tpp = newp;
  22. ts = 1; // height of toroid
  23. tx = radius; // radial width of toroid, measured in centre of ring
  24. tl = 0; // centre of rotation
  25. Point(tpp ) = { tx, 0, 0, lc};
  26. Point(tpp+1) = { ts+tx, 0, 0, lc};
  27. Point(tpp+2) = { tx, ts, 0, lc};
  28. Point(tpp+3) = { tx, -ts, 0, lc};
  29. Point(tpp+4) = {-ts+tx, 0, 0, lc};
  30. cc = newc;
  31. Circle(cc ) = {tpp+1, tpp, tpp+2};
  32. Circle(cc+1) = {tpp+2, tpp, tpp+4};
  33. Circle(cc+2) = {tpp+4, tpp, tpp+3};
  34. Circle(cc+3) = {tpp+3, tpp, tpp+1};
  35. ll = newll;
  36. Line Loop(ll) = {cc, cc+1, cc+2, cc+3};
  37. ps = news;
  38. pio2=Pi/2;
  39. Plane Surface(ps) = {ll};
  40. tv1[] = Extrude {{0, 1, 0},{-tl,0,0}, 2*Pi/3} { Surface{ps}; };
  41. tv2[] = Extrude {{0, 1, 0},{-tl,0,0}, 2*Pi/3} { Surface{28}; };
  42. tv3[] = Extrude {{0, 1, 0},{-tl,0,0}, 2*Pi/3} { Surface{50}; };
  43. //t1[] = Rotate {{0,0,1},{0,0,0},pio2 } {Duplicata{Surface{ps};}};
  44. //Extrude Surface {ps, {0,1,0}, {-tl,0,0}, 2*Pi/3} { Recombine ;};
  45. //Extrude Surface {28, {0,1,0}, {-tl,0,0}, 2*Pi/3}; //{Layers{10,73,1};};
  46. //Extrude Surface {50, {0,1,0}, {-tl,0,0}, 2*Pi/3}; //{Layers{10,73,1};};
  47. /* Make a list of a ring (annulus) of surfaces around the hole */
  48. allParts[] = {tv1[0], tv2[0], tv3[0]};
  49. /* Make surfaces to be meshed by transfinite algorithm */
  50. //Transfinite Surface {allParts[]};
  51. /* The "Recombine Surface" command is issued in order to
  52. * crate quadrilateral elements.
  53. */
  54. //Recombine Surface {allParts[]};
  55. // Extrude Surface {12, {0,0,1}, {0,0,0}, 2*Pi/3} {
  56. // Recombine ; Layers { 6, 54, 1 } ;
  57. // } ;
  58. // Total Solution Space
  59. X0 = -Box;
  60. X1 = Box;
  61. Y0 = -Box;
  62. Y1 = Box;
  63. Z0 = -Box;
  64. Z1 = Box;
  65. /////////////////////////////////////
  66. // Large Bounding box
  67. pp = newp;
  68. Point(pp) = {X0, Y0, Z0, blc};
  69. Point(pp+1) = {X1, Y0, Z0, blc};
  70. Point(pp+2) = {X1, Y1, Z0, blc};
  71. Point(pp+3) = {X0, Y1, Z0, blc};
  72. //
  73. lv = newl;
  74. Line(lv) = {pp,pp+1};
  75. Line(lv+1) = {pp+1,pp+2};
  76. Line(lv+2) = {pp+2,pp+3};
  77. Line(lv+3) = {pp+3,pp};
  78. Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
  79. //
  80. // Hard coded doom
  81. bs = news;
  82. Plane Surface(bs) = {lv+4};
  83. //
  84. //v = newv;
  85. v[] = Extrude {0, 0, Z1-Z0} { Surface{bs}; };
  86. /* This is GOOD */
  87. Surface{ allParts } In Volume{v[1]};
  88. //Surface{t1[0]} In Volume{v[1]};
  89. //Surface{t2[0]} In Volume{v[1]};
  90. //Surface{t3[0]} In Volume{v[1]};