Galerkin FEM for elliptic PDEs
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120
  1. /**
  2. * Trevor P. Irons, XRI Geophysics, LLC
  3. * Test of Coloumbic Magnetic Potential
  4. */
  5. lc = 1e-2; // Target element size
  6. R = .05; // Magnet Radius
  7. D0 = 9.5; // Top of mesh
  8. D1 = 11.5; // Bottom of mesh
  9. // Total Solution Space
  10. X0 = -2.;
  11. X1 = 2.;
  12. Y0 = -2.;
  13. Y1 = 2.;
  14. //Z0 = 9.5;
  15. //Z1 = 11.5;
  16. ////////////////////////////////////
  17. // North Pole
  18. Function Ring
  19. X *= 1.2;
  20. // North Pole
  21. p0 = newp; Point(p0) = { 0, 0, D0, lc};
  22. p1 = newp; Point(p1) = { X*R, 0, D0, lc};
  23. p2 = newp; Point(p2) = { 0, X*R, D0, lc};
  24. p3 = newp; Point(p3) = {-X*R, 0, D0, lc};
  25. c1 = newc; Circle(c1) = {p1, p0, p3};
  26. l1 = newl; Line(l1) = {p3,p1};
  27. l2 = newl; Line Loop(l2) = {c1, l1};
  28. s1 = news; Plane Surface(s1) = {l2};
  29. Extrude {0, 0, D1-D0} { Surface{s1}; }
  30. // South Pole
  31. p4 = newp; Point(p4) = { 0, -0.0001, D0, lc}; // 8
  32. p5 = newp; Point(p5) = { X*R, -0.0001, D0, lc}; // 9
  33. p6 = newp; Point(p6) = { 0, -X*R, D0, lc}; // 10
  34. p7 = newp; Point(p7) = {-X*R, -0.0001, D0, lc}; // 11
  35. c2 = newc; Circle(c2) = {p7, p4, p5};
  36. l3 = newl; Line(l3) = {p5,p7};
  37. l4 = newl; Line Loop(l4) = {c2, l3};
  38. s2 = news; Plane Surface(s2) = {l4};
  39. Extrude {0, 0, D1-D0} { Surface{s2}; }
  40. //Compound Volume(4) = {vol2[1], vol3[1]};
  41. //Circle(12) = {11, 8, 9};
  42. //Line(13) = {9, 11};
  43. //Line Loop(14) = {12, 13};
  44. //Plane Surface(15) = {14};
  45. Return
  46. X = 1.;
  47. For r In {1:10}
  48. Call Ring;
  49. EndFor
  50. /*
  51. Point(0) = { 0, 0, D0, lc};
  52. Point(1) = { R, 0, D0, lc};
  53. Point(2) = { 0, R, D0, lc};
  54. Point(3) = {-R, 0, D0, lc};
  55. // Connect up the points
  56. Circle(4) = {1, 0, 3};
  57. Line(5) = {3,1};
  58. Line Loop(6) = {4,5};
  59. Plane Surface(7) = {6};
  60. ////////////////////////////////////
  61. // South Pole
  62. Point(8) = { 0, -0.0001, D0, lc};
  63. Point(9) = { R, -0.0001, D0, lc};
  64. Point(10) = { 0, -R, D0, lc};
  65. Point(11) = {-R, -0.0001, D0, lc};
  66. // Connect up the points
  67. Circle(12) = {11, 8, 9};
  68. Line(13) = {9, 11};
  69. Line Loop(14) = {12, 13};
  70. Plane Surface(15) = {14};
  71. //////////////////////////////////////
  72. // Extrude magnet
  73. Extrude {0, 0, D1-D0} { Surface{7}; }
  74. Extrude {0, 0, D1-D0} { Surface{15}; }
  75. /////////////////////////////////////
  76. // Large Bounding box
  77. Point(116) = {X0, Y0, Z0};
  78. Point(117) = {X1, Y0, Z0};
  79. Point(118) = {X1, Y1, Z0};
  80. Point(119) = {X0, Y1, Z0};
  81. Line(120) = {116,117};
  82. Line(121) = {117,118};
  83. Line(122) = {118,119};
  84. Line(123) = {119,116};
  85. Line Loop(124) = {120, 121, 122, 123};
  86. Plane Surface(125) = {124};
  87. Extrude {0, 0, Z1-Z0} { Surface{125}; }
  88. //////////////////////////////////////
  89. // Rings of sensitivity calculation
  90. //Circle(222) = {2.*1, 0, 2.*2};
  91. //////////////////////////////////////
  92. // Volumes
  93. Physical Volume("plus") = {2};
  94. Physical Volume("minus") = {1};
  95. Physical Volume("background") = {3};
  96. Coherence;
  97. */