Galerkin FEM for elliptic PDEs
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198
  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 10/31/2014 10:44:16 AM
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@xri-geo.com
  14. * @copyright Copyright (c) 2014, XRI Geophysics, LLC
  15. * @copyright Copyright (c) 2014, Trevor Irons
  16. */
  17. /*********************************************************************
  18. * Starting model was
  19. * Gmsh tutorial 4, modified to give tunnel model with simple topography
  20. *
  21. * Built-in functions, holes, strings, mesh color
  22. *
  23. *********************************************************************/
  24. // As usual, we start by defining some variables:
  25. _m = 1; // _m characteristic length thingy change to m if probably
  26. // Distance
  27. e1 = 50. * _m; // Northing length of the whole model (twice this number)
  28. e2 = 5. * _m ; // Width of first dip
  29. h1 = 15.5 * _m; // height of Low wall
  30. h2 = 50 * _m; // Height of the Low wall
  31. h5 = 2 * _m; // Height of the Dip
  32. h3 = -8 * _m; // Depth of tunnel
  33. h4 = 2 * _m; // Height of tunnel
  34. he = 100*_m; // Extrusion lengths
  35. sy = -he/2.;
  36. R1 = 1 * _m; // Dip, extremely fragile!
  37. R2 = 1. * _m; // Tunnel width
  38. // Define Electrode grid
  39. NE = 2; // Num easting
  40. NN = 40; // Num northing
  41. SE = -10;
  42. SN = -20;
  43. de = 1;
  44. dn = 1;
  45. Lc1 = 5.5 * _m; // Characteristic length of thingies
  46. Lc2 = 5.5 * _m;
  47. //Lc3 = 5.5 * _m;
  48. Lc4 = 0.15*de * _m; // Electrode sensitivity, 1/4 electrode spacing is good
  49. //Lc4 = 5.5*de * _m; // Electrode sensitivity, 1/4 electrode spacing is good
  50. // We can use all the usual mathematical functions (note the
  51. // capitalized first letters), plus some useful functions like
  52. // Hypot(a, b) := Sqrt(a^2 + b^2):
  53. ccos = (-h5*R1 + e2 * Hypot(h5, Hypot(e2, R1))) / (h5^2 + e2^2);
  54. ssin = Sqrt(1 - ccos^2);
  55. p1 = newp;
  56. // Then we define some points and some lines using these variables:
  57. Point(p1 ) = {-e1-e2, sy, 0 , Lc1};
  58. Point(p1+1) = {-e1-e2, sy, h1+h2, Lc1};
  59. Point(p1+2) = {e1+e2 , sy, h1+h2, Lc1};
  60. Point(p1+3) = { e1+e2, sy, h1 , Lc1};
  61. // Defines the dip
  62. Point(p1+4)= { R1 / ssin, sy, h5+R1*ccos, Lc2}; // p5
  63. Point(p1+5)= {-R1 / ssin, sy, h5+R1*ccos, Lc2}; // p6
  64. // Defines something else
  65. Point(p1+6) = {-e2 , sy, 0.0, Lc2};
  66. /* Tunnel */ /*
  67. Point(p1+7) = {-R2 , sy, h1+h3 , Lc2};
  68. Point(p1+8) = {-R2 , sy, h1+h3+h4, Lc2};
  69. Point(p1+9) = { 0 , sy, h1+h3+h4, Lc2};
  70. Point(p1+10) = { R2 , sy, h1+h3+h4, Lc2};
  71. Point(p1+11) = { R2 , sy, h1+h3 , Lc2};
  72. Point(p1+12) = { 0 , sy, h1+h3 , Lc2};
  73. Point(p1+13) = { 0 , sy, h1+h3+h4+R2, Lc2};
  74. Point(p1+14) = { 0 , sy, h1+h3-R2, Lc2};
  75. */
  76. L = newl;
  77. Line(L ) = {p1 , p1+6};
  78. Line(L+1) = {p1+6, p1+5};
  79. Line(L+2) = {p1+5, p1+4};
  80. Line(L+3) = {p1+4, p1+3};
  81. Line(L+4) = {p1+3, p1+2}; // MAYBE SHOULD BE 11,6??
  82. Line(L+5) = {p1+2, p1+1};
  83. Line(L+6) = {p1+1, p1 };
  84. /*
  85. Line(L+7) = {p1+8, p1+7}; // TUNNEL
  86. Line(L+8) = {p1+11, p1+10};
  87. C = newc;
  88. Circle(C ) = {p1+10, p1+9, p1+13};
  89. Circle(C+1) = {p1+13, p1+9, p1+8 };
  90. Circle(C+2) = {p1+7, p1+12, p1+14};
  91. Circle(C+3) = {p1+14, p1+12, p1+11};
  92. */
  93. LL = newll;
  94. //Line Loop(LL) = {C+1, L+7, C+2, C+3, L+8, C}; // THIS IS THE TUNNEL
  95. Line Loop(LL+1) = { L, L+1, L+2, L+3, L+4, L+5, L+6 };
  96. s1 = news;
  97. //Plane Surface(s1) = {LL+1, LL};
  98. Plane Surface(s1) = {LL+1};
  99. // Extrude
  100. out[] = Extrude {0,he,0} {
  101. Surface{s1};
  102. };
  103. //////////////////////////////////////////////////////////
  104. // Dirichlet Boundary
  105. // We don't want these to be physical surfaces, BUT, we can use them to embed point into...
  106. // Current workflow requires meshing twice, once with this and then saving .stl file then again as a vtk without.
  107. Physical Surface(2) = {s1, out[0], out[6], out[7], out[8] }; //
  108. // Mesh Volume
  109. Physical Volume(1) = {out[1]};
  110. ep = newp;
  111. Printf("%YAML 1.2") > "electrodes.yaml";
  112. Printf("---") >> "electrodes.yaml";
  113. Printf("Electrodes:") >> "electrodes.yaml";
  114. ii = 0;
  115. For iie In {0:NE-1}
  116. For iin In {0:NN-1}
  117. yy = SE+de*iie; // easting
  118. xx = SN+dn*iin; // northing
  119. zz = 0;
  120. If (xx < -e2) // Point on top shelf
  121. Point(ep+ii) = {xx , yy, 0.0, Lc4};
  122. Point{ep+ii} In Surface{ out[2] }; // out[2] is top shelf
  123. EndIf
  124. If ( xx>-e2 && xx < -R1 / ssin )
  125. zz = (xx+e2) * ( (h5+R1*ccos) / ((e2) - (R1/ssin)) );
  126. Point(ep+ii) = {xx, yy, zz, Lc4};
  127. Point{ep+ii} In Surface{ out[3] }; // out[3] is next shelf
  128. EndIf
  129. If( xx> -R1 / ssin && xx< R1/ssin)
  130. zz = h5+R1*ccos;
  131. Point(ep+ii) = {xx , yy, zz, Lc4};
  132. Point{ep+ii} In Surface{ out[4] }; // out[4] is next shelf
  133. EndIf
  134. If( xx>R1/ssin )
  135. zz = h5+R1*ccos + ( (xx-(R1/ssin)) * ( (h1-(h5+R1*ccos)) / ((e1+e2) - (R1/ssin)) ));
  136. Point(ep+ii) = {xx , yy, zz, Lc4};
  137. Point{ep+ii} In Surface{ out[5] }; // out[5] is next shelf
  138. EndIf
  139. // TODO 26 is a magic number for this file. Seems to work but check
  140. Printf(" L%g-%g: !<DCIPElectrode> &L%g-%g", iie, iin, iie, iin) >> "electrodes.yaml";
  141. Printf(" Node_ID: %g", 26+ii) >> "electrodes.yaml";
  142. Printf(" Location: !<Vector3r>") >> "electrodes.yaml";
  143. Printf(" -") >> "electrodes.yaml";
  144. Printf(" - %f", xx) >> "electrodes.yaml";
  145. Printf(" - %f", yy) >> "electrodes.yaml";
  146. Printf(" - %f", zz) >> "electrodes.yaml";
  147. //Physical Point(ii) = {ep+ii};
  148. ii += 1;
  149. EndFor
  150. EndFor
  151. ///////////////////////////////////////////////
  152. // Attractor Field
  153. Field[1] = Attractor;
  154. Field[1].NodesList = {ep:ep+NE*NN}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
  155. Field[2] = Threshold;
  156. Field[2].IField = 1;
  157. Field[2].LcMin = Lc4;
  158. Field[2].LcMax = Lc1;
  159. Field[2].DistMin = .1;
  160. Field[2].DistMax = 20;
  161. // Use minimum of all the fields as the background field
  162. Field[3] = Min;
  163. Field[3].FieldsList = {2};
  164. Background Field = 3;