FEM4EllipticPDE/examples/DC/tunnel.geo

/* 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: !<DCIPElectrode> &L%g-%g", iie, iin, iie, iin) >> "electrodes.yaml";
      Printf("    Node_ID: %g", 26+ii) >> "electrodes.yaml";
      Printf("    Location: !<Vector3r>") >> "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;