/* This file is part of Lemma, a geophysical modelling and inversion API. * More information is available at http://lemmasoftware.org */ /* This Source Code Form is subject to the terms of the Mozilla Public * License, v. 2.0. If a copy of the MPL was not distributed with this * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ /** * @file * @date 10/31/2014 10:44:16 AM * @version $Id$ * @author Trevor Irons (ti) * @email Trevor.Irons@xri-geo.com * @copyright Copyright (c) 2014, XRI Geophysics, LLC * @copyright Copyright (c) 2014, Trevor Irons */ /********************************************************************* * Starting model was * Gmsh tutorial 4, modified to give tunnel model with simple topography * * Built-in functions, holes, strings, mesh color * *********************************************************************/ // As usual, we start by defining some variables: _m = 1; // _m characteristic length thingy change to m if probably // Distance e1 = 50. * _m; // Northing length of the whole model (twice this number) e2 = 5. * _m ; // Width of first dip h1 = 15.5 * _m; // height of Low wall h2 = 50 * _m; // Height of the Low wall h5 = 2 * _m; // Height of the Dip h3 = -8 * _m; // Depth of tunnel h4 = 2 * _m; // Height of tunnel he = 100*_m; // Extrusion lengths sy = -he/2.; R1 = 1 * _m; // Dip, extremely fragile! R2 = 1. * _m; // Tunnel width // Define Electrode grid NE = 2; // Num easting NN = 40; // Num northing SE = -10; SN = -20; de = 1; dn = 1; Lc1 = 5.5 * _m; // Characteristic length of thingies Lc2 = 5.5 * _m; //Lc3 = 5.5 * _m; Lc4 = 0.15*de * _m; // Electrode sensitivity, 1/4 electrode spacing is good //Lc4 = 5.5*de * _m; // Electrode sensitivity, 1/4 electrode spacing is good // We can use all the usual mathematical functions (note the // capitalized first letters), plus some useful functions like // Hypot(a, b) := Sqrt(a^2 + b^2): ccos = (-h5*R1 + e2 * Hypot(h5, Hypot(e2, R1))) / (h5^2 + e2^2); ssin = Sqrt(1 - ccos^2); p1 = newp; // Then we define some points and some lines using these variables: Point(p1 ) = {-e1-e2, sy, 0 , Lc1}; Point(p1+1) = {-e1-e2, sy, h1+h2, Lc1}; Point(p1+2) = {e1+e2 , sy, h1+h2, Lc1}; Point(p1+3) = { e1+e2, sy, h1 , Lc1}; // Defines the dip Point(p1+4)= { R1 / ssin, sy, h5+R1*ccos, Lc2}; // p5 Point(p1+5)= {-R1 / ssin, sy, h5+R1*ccos, Lc2}; // p6 // Defines something else Point(p1+6) = {-e2 , sy, 0.0, Lc2}; /* Tunnel */ /* Point(p1+7) = {-R2 , sy, h1+h3 , Lc2}; Point(p1+8) = {-R2 , sy, h1+h3+h4, Lc2}; Point(p1+9) = { 0 , sy, h1+h3+h4, Lc2}; Point(p1+10) = { R2 , sy, h1+h3+h4, Lc2}; Point(p1+11) = { R2 , sy, h1+h3 , Lc2}; Point(p1+12) = { 0 , sy, h1+h3 , Lc2}; Point(p1+13) = { 0 , sy, h1+h3+h4+R2, Lc2}; Point(p1+14) = { 0 , sy, h1+h3-R2, Lc2}; */ L = newl; Line(L ) = {p1 , p1+6}; Line(L+1) = {p1+6, p1+5}; Line(L+2) = {p1+5, p1+4}; Line(L+3) = {p1+4, p1+3}; Line(L+4) = {p1+3, p1+2}; // MAYBE SHOULD BE 11,6?? Line(L+5) = {p1+2, p1+1}; Line(L+6) = {p1+1, p1 }; /* Line(L+7) = {p1+8, p1+7}; // TUNNEL Line(L+8) = {p1+11, p1+10}; C = newc; Circle(C ) = {p1+10, p1+9, p1+13}; Circle(C+1) = {p1+13, p1+9, p1+8 }; Circle(C+2) = {p1+7, p1+12, p1+14}; Circle(C+3) = {p1+14, p1+12, p1+11}; */ LL = newll; //Line Loop(LL) = {C+1, L+7, C+2, C+3, L+8, C}; // THIS IS THE TUNNEL Line Loop(LL+1) = { L, L+1, L+2, L+3, L+4, L+5, L+6 }; s1 = news; //Plane Surface(s1) = {LL+1, LL}; Plane Surface(s1) = {LL+1}; // Extrude out[] = Extrude {0,he,0} { Surface{s1}; }; ////////////////////////////////////////////////////////// // Dirichlet Boundary // We don't want these to be physical surfaces, BUT, we can use them to embed point into... // Current workflow requires meshing twice, once with this and then saving .stl file then again as a vtk without. Physical Surface(2) = {s1, out[0], out[6], out[7], out[8] }; // // Mesh Volume Physical Volume(1) = {out[1]}; ep = newp; Printf("%YAML 1.2") > "electrodes.yaml"; Printf("---") >> "electrodes.yaml"; Printf("Electrodes:") >> "electrodes.yaml"; ii = 0; For iie In {0:NE-1} For iin In {0:NN-1} yy = SE+de*iie; // easting xx = SN+dn*iin; // northing zz = 0; If (xx < -e2) // Point on top shelf Point(ep+ii) = {xx , yy, 0.0, Lc4}; Point{ep+ii} In Surface{ out[2] }; // out[2] is top shelf EndIf If ( xx>-e2 && xx < -R1 / ssin ) zz = (xx+e2) * ( (h5+R1*ccos) / ((e2) - (R1/ssin)) ); Point(ep+ii) = {xx, yy, zz, Lc4}; Point{ep+ii} In Surface{ out[3] }; // out[3] is next shelf EndIf If( xx> -R1 / ssin && xx< R1/ssin) zz = h5+R1*ccos; Point(ep+ii) = {xx , yy, zz, Lc4}; Point{ep+ii} In Surface{ out[4] }; // out[4] is next shelf EndIf If( xx>R1/ssin ) zz = h5+R1*ccos + ( (xx-(R1/ssin)) * ( (h1-(h5+R1*ccos)) / ((e1+e2) - (R1/ssin)) )); Point(ep+ii) = {xx , yy, zz, Lc4}; Point{ep+ii} In Surface{ out[5] }; // out[5] is next shelf EndIf // TODO 26 is a magic number for this file. Seems to work but check Printf(" L%g-%g: ! &L%g-%g", iie, iin, iie, iin) >> "electrodes.yaml"; Printf(" Node_ID: %g", 26+ii) >> "electrodes.yaml"; Printf(" Location: !") >> "electrodes.yaml"; Printf(" -") >> "electrodes.yaml"; Printf(" - %f", xx) >> "electrodes.yaml"; Printf(" - %f", yy) >> "electrodes.yaml"; Printf(" - %f", zz) >> "electrodes.yaml"; //Physical Point(ii) = {ep+ii}; ii += 1; EndFor EndFor /////////////////////////////////////////////// // Attractor Field Field[1] = Attractor; 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}; Field[2] = Threshold; Field[2].IField = 1; Field[2].LcMin = Lc4; Field[2].LcMax = Lc1; Field[2].DistMin = .1; Field[2].DistMax = 20; // Use minimum of all the fields as the background field Field[3] = Min; Field[3].FieldsList = {2}; Background Field = 3;