Galerkin FEM for elliptic PDEs
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

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]};