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.

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485
  1. /**
  2. * Trevor P. Irons, XRI Geophysics, LLC
  3. * Test of Coloumbic Magnetic Potential
  4. */
  5. lc = 1e-2; // Target element size
  6. R = .11; // Minimum Radius
  7. R2 = .175; // Minimum Radius
  8. D0 = 9.75; // Top of mesh
  9. D1 = 11.25; // Bottom of mesh
  10. ////////////////////////////////////
  11. // North Pole
  12. Function Ring
  13. // centre
  14. p0 = newp; Point(p0) = { 0, 0, D0, lc};
  15. // Points defining outer ring
  16. p1 = newp; Point(p1) = { R, 0, D0, lc};
  17. p2 = newp; Point(p2) = { 0, R, D0, lc};
  18. p3 = newp; Point(p3) = {-R, 0, D0, lc};
  19. p4 = newp; Point(p4) = { 0, -R, D0, lc};
  20. c1 = newc; Circle(c1) = {p1, p0, p2};
  21. c2 = newc; Circle(c2) = {p2, p0, p3};
  22. c3 = newc; Circle(c3) = {p3, p0, p4};
  23. c4 = newc; Circle(c4) = {p4, p0, p1};
  24. // Inner ring
  25. p5 = newp; Point(p5) = { R2, 0, D0, lc};
  26. p6 = newp; Point(p6) = { 0, R2, D0, lc};
  27. p7 = newp; Point(p7) = {-R2, 0, D0, lc};
  28. p8 = newp; Point(p8) = { 0, -R2, D0, lc};
  29. c5 = newc; Circle(c5) = {p5, p0, p6};
  30. c6 = newc; Circle(c6) = {p6, p0, p7};
  31. c7 = newc; Circle(c7) = {p7, p0, p8};
  32. c8 = newc; Circle(c8) = {p8, p0, p5};
  33. l1 = newl; Line Loop(l1) = {c1, c2, c3, c4, c5, c6, c7, c8};
  34. s1 = news; Plane Surface(s1) = {l1};
  35. Extrude {0,0,D1-D0} {
  36. Surface{s1};
  37. }
  38. Return
  39. Function Cyl
  40. // centre
  41. p0 = newp; Point(p0) = { 0, 0, D0, lc};
  42. // Points defining outer ring
  43. p1 = newp; Point(p1) = { R2, 0, D0, lc};
  44. p2 = newp; Point(p2) = { 0, R2, D0, lc};
  45. p3 = newp; Point(p3) = {-R2, 0, D0, lc};
  46. p4 = newp; Point(p4) = { 0, -R2, D0, lc};
  47. c1 = newc; Circle(c1) = {p1, p0, p2};
  48. c2 = newc; Circle(c2) = {p2, p0, p3};
  49. c3 = newc; Circle(c3) = {p3, p0, p4};
  50. c4 = newc; Circle(c4) = {p4, p0, p1};
  51. l1 = newl; Line Loop(l1) = {c1, c2, c3, c4};
  52. s1 = news; Plane Surface(s1) = {l1};
  53. Extrude {0,0,D1-D0} {
  54. Surface{s1};
  55. }
  56. Return
  57. Call Cyl;
  58. /*
  59. For r In {1:12}
  60. Call Ring;
  61. R2 = R;
  62. R *= 1.2;
  63. lc *= 1.1;
  64. EndFor
  65. */