Examples now build for FEM and install to installdir/share

CMakeLists.txt

@@ -5,3 +5,9 @@ add_subdirectory("src")
 
 add_library( fem4ellipticpde ${FEMSOURCE} )  
 target_link_libraries(fem4ellipticpde "lemmacore")
+
+install ( TARGETS fem4ellipticpde DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
+
+if (LEMMABUILDEXAMPLES)
+	add_subdirectory(examples)
+endif()

examples/CMakeLists.txt

@@ -0,0 +1,33 @@
+add_executable( FEM4EllipticPDE_bhmag FEM4EllipticPDE_bhmag.cpp  )
+target_link_libraries(  FEM4EllipticPDE_bhmag  "lemmacore" "fem4ellipticpde")
+
+add_executable( FEM4EllipticPDE FEM4EllipticPDE.cpp  )
+target_link_libraries(  FEM4EllipticPDE  "lemmacore" "fem4ellipticpde")
+
+add_executable( merge merge.cpp  )
+target_link_libraries(  merge  "lemmacore" "fem4ellipticpde")
+
+add_executable( VTKDC VTKDC.cpp  )
+target_link_libraries( VTKDC  "lemmacore" "fem4ellipticpde")
+
+add_executable( VTKEdgeG VTKEdgeG.cpp  )
+target_link_libraries( VTKEdgeG  "lemmacore" "fem4ellipticpde")
+
+add_executable( VTKEdgeGsphere VTKEdgeGsphere.cpp  )
+target_link_libraries( VTKEdgeGsphere  "lemmacore" "fem4ellipticpde")
+
+#INSTALL_TARGETS( "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
+INSTALL_TARGETS( "/share/FEM4EllipticPDE/"
+	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKEdgeGsphere
+)
+
+install( DIRECTORY "borehole"
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
+	PATTERN .svn EXCLUDE
+)
+
+install( DIRECTORY "DC"
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
+	PATTERN .svn EXCLUDE
+)
+

examples/DC/DC.yaml

@@ -0,0 +1,126 @@
+!<DCSurvey>
+NumberOfElectrodes: 9
+Injections:
+  J-0:
+    A: !<DCIPElectrode> &1
+      Node_ID: -1
+      Label: L0-E0
+      Location: !<Vector3r>
+        -
+          - 13.23
+          - 18.1
+          - -25
+    B: !<DCIPElectrode> &2
+      Location: !<Vector3r>
+        -
+          - 13.23
+          - 18.1
+          - -25
+      Node_ID: -1
+      Label: L0-E1
+    J: 1
+    Measurements:
+      Number: 4
+      V-0:
+        M: !<DCIPElectrode> &3
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E2
+        N: !<DCIPElectrode> &4
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E3
+      V-1:
+        M: !<DCIPElectrode> &5
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E4
+        N: !<DCIPElectrode> &6
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E5
+      V-2:
+        M: *6
+        N: !<DCIPElectrode> &7
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E6
+      V-3:
+        M: !<DCIPElectrode> &8
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E7
+        N: !<DCIPElectrode> &9
+          Location: !<Vector3r>
+            -
+              - 13.23
+              - 18.1
+              - -25
+          Node_ID: -1
+          Label: L0-E8
+  J-1:
+    A: *2
+    B: *3
+    J: 1
+    Measurements:
+      Number: 3
+      V-0:
+        M: *3
+        N: *4
+      V-1:
+        M: *5
+        N: *6
+      V-2:
+        M: *6
+        N: *7
+  J-2:
+    A: *3
+    B: *4
+    J: 1
+    Measurements:
+      Number: 3
+      V-0:
+        M: *3
+        N: *4
+      V-1:
+        M: *5
+        N: *6
+      V-2:
+        M: *6
+        N: *7
+Electrodes:
+  L0-E2: *3
+  L0-E3: *4
+  L0-E4: *5
+  L0-E5: *6
+  L0-E6: *7
+  L0-E7: *8
+  L0-E8: *9
+  L0-E0: *1
+  L0-E1: *2
+---
+

examples/DC/tunnel.geo

@@ -0,0 +1,198 @@
+/* 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;
+

examples/borehole/Makefile

@@ -0,0 +1,27 @@
+
+sphereG.vtu : sphere.geo
+	#gmsh  -3 -format vtk -o sphere.vtk  sphere.geo
+	../utVTKEdgeGsphere sphere.vtk .25 gravity  sphereG.vtu
+	#../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
+
+sphereOut: 
+	../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 
+	#paraview --state=sliceSphere.pvsm
+
+all : sphereG.vtu sphereOut
+
+# --data=sphereOut.vtu
+
+# even finer meshing
+# radius = 4, u = \pm .9              analytic = \pm 1.33
+
+# Finer meshing
+# radius = 5, u = \pm 1.65            analytic = \pm 1.67
+# radius = 4, u = \pm 1.09            analytic = \pm 1.33
+
+# radius = 5, u = \pm 2.5             analytic = \pm 1.67
+# radius = 4, u = \pm 1.66            analytic = \pm 1.33
+# radius = 3, u= \pm .926 -.985       analytic = \pm 1
+# radius = 2, u = \pm.42              analytic = \pm .667
+# for radius=1,  u = \pm .111         analytic = \pm .333 
+# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167

examples/borehole/Notes.txt

@@ -0,0 +1,31 @@
+Gmsh does not seem to import the group into exported VTK dataset. Therefore the following workflow seems to be 
+needed.
+
+In Gmsh, make a seperate .geo file for each body. Then export the grid for each physical volume into VTK format.
+Then you load all of these in Paraview, Use Calculator to assign values, and then use AppendDatasets. Then you 
+can save as a .vtu file. 
+
+For just the grid, you can save the whole Gmsh grid as a VTK grid. 
+
+=================================================================================================================
+
+Since we only want \nabla \cdot M on the calculation mesh do the following:
+
+	I.
+		1: In Gmsh make .vtk files for Grid_merge 
+		2: In Gmsh make .vtk files for BoreholeMagnetPlus (2D grid!) 
+			When doing this make sure you save all in the GUI. Or else no cells are exported.     
+ 
+	II. In Paraview
+		1. Open Grid_merge.vtk file
+		2. Open BoreholeMagnetPlus.vtk file
+			->Use Calculator assign 0 to everything
+			->Clip with normal
+				Use Calculator to assign -1
+			-> Clip with -normal 
+				Use Calculator to assign 1
+			-> Clip of ends. 
+		3. cntrl-click the three Calculator datasets. And use AppendDatasets Filter.
+		4. Save Data (Select Appended) as .vtu
+		5. Disco  		
+

examples/borehole/boreholeMagnet.geo

@@ -0,0 +1,100 @@
+/**
+ *  Trevor P. Irons, XRI Geophysics, LLC
+ *  Test of Coloumbic Magnetic Potential
+ */
+
+lc =   .01;   // Target element size
+R =    .25;   // Magnet Radius
+D0 = 10;      // Top of magnet
+D1 = 11;      // Bottom of magnet
+
+// Total Solution Space
+X0 = -1.;
+X1 =  1.;
+Y0 = -1.;
+Y1 =  1.;
+Z0 =  9.5;
+Z1 = 11.5;
+
+////////////////////////////////////
+// North Pole
+p = newp;
+dd = 5;
+d = lc/(1*dd);
+Point(p) =   { 0,  0, D0, lc/dd};
+Point(p+1) = { R,  d, D0, lc/dd};
+Point(p+2) = { R, -d, D0, lc/dd};
+Point(p+3) = {-R,  d, D0, lc/dd};
+Point(p+4) = {-R, -d, D0, lc/dd};
+
+// Connect up the points
+c = newc;
+Circle(c)  = {p+3, p, p+1};
+
+c2 = newc;
+Circle(c2)  = {p+4, p, p+2};
+
+ss = news;
+ss = Extrude {0, 0, D1-D0} { Line{c}; };
+
+ss2 = news;
+ss2 = Extrude {0, 0, D1-D0} { Line{c2}; };
+
+
+//Symmetry { 1.0, 1.0, 1.0, 0.0 }{Duplicata{Surface{ss[1]};}}
+//lcs = newl;
+//lcs = Symmetry { 0.0, 0.0,.0001, 0.0 }{Duplicata{Line{c};}};
+//lcs =
+//Symmetry { 0.0, 1.0, 0.0, 0.0 }{Duplicata{Line{c};}}
+
+//Symmetry {0, 0, 120, 0} {
+// Duplicata{Surface{ss[1]};}
+//}
+
+
+
+
+//ss2 = news;
+//ss = Extrude {0, 0, D1-D0} { Line{c, lcs}; };
+//ss2 = Extrude {0, 0, D1-D0} { Line{l, lcs}; };
+
+//Symmetry {1,0,0,0}{ Duplicata{Surface{ss[1]};}};
+//Symmetry{ expression-list }{ transform-list }:
+
+//ss2 = news;
+//ss2 = Extrude {0, 0, D1-D0} { Line{c2, l2}; };
+
+/////////////////////////////////////
+// Large Bounding box
+pp = newp;
+Point(pp)    = {X0, Y0, Z0, lc};
+Point(pp+1)  = {X1, Y0, Z0, lc};
+Point(pp+2)  = {X1, Y1, Z0, lc};
+Point(pp+3)  = {X0, Y1, Z0, lc};
+
+
+lv = newl;
+Line(lv) = {pp,pp+1};
+Line(lv+1) = {pp+1,pp+2};
+Line(lv+2) = {pp+2,pp+3};
+Line(lv+3) = {pp+3,pp};
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
+
+// Hard coded doom
+Plane Surface(125) = {lv+4};
+
+//v = newv;
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
+
+Surface{ss[1]} In Volume{v[1]};
+Surface{ss2[1]} In Volume{v[1]};
+
+Physical Volume(1) = {v[1]};
+
+//Physical Volume("minus") = {1};
+//Physical Volume("background") = {3};
+//Coherence;
+
+//Coherence;
+
+//Mesh.RecombineAll = 1;

examples/borehole/boreholeMagnetGrid.geo

@@ -0,0 +1,120 @@
+/**
+ *  Trevor P. Irons, XRI Geophysics, LLC
+ *  Test of Coloumbic Magnetic Potential
+ */
+
+lc = 1e-2;   // Target element size
+R =  .05;    // Magnet Radius
+D0 = 9.5;    // Top of mesh
+D1 = 11.5;   // Bottom of mesh
+
+// Total Solution Space
+X0 = -2.;
+X1 =  2.;
+Y0 = -2.;
+Y1 =  2.;
+//Z0 =  9.5;
+//Z1 = 11.5;
+
+////////////////////////////////////
+// North Pole
+
+Function Ring
+    X *= 1.2;
+    // North Pole
+    p0 = newp; Point(p0) = {   0,   0, D0, lc};
+    p1 = newp; Point(p1) = { X*R,   0, D0, lc};
+    p2 = newp; Point(p2) = {   0, X*R, D0, lc};
+    p3 = newp; Point(p3) = {-X*R,   0, D0, lc};
+
+    c1 = newc; Circle(c1) = {p1, p0, p3};
+    l1 = newl; Line(l1) = {p3,p1};
+    l2 = newl; Line Loop(l2) = {c1, l1};
+    s1 = news; Plane Surface(s1) = {l2};
+    Extrude {0, 0, D1-D0} { Surface{s1};  }
+
+    // South Pole
+    p4 = newp; Point(p4)  = {   0,  -0.0001, D0, lc}; //  8
+    p5 = newp; Point(p5)  = { X*R,  -0.0001, D0, lc}; //  9
+    p6 = newp; Point(p6)  = {   0,  -X*R,    D0, lc}; // 10
+    p7 = newp; Point(p7)  = {-X*R,  -0.0001, D0, lc}; // 11
+
+    c2 = newc; Circle(c2) = {p7, p4, p5};
+    l3 = newl; Line(l3) = {p5,p7};
+    l4 = newl; Line Loop(l4) = {c2, l3};
+    s2 = news; Plane Surface(s2) = {l4};
+    Extrude {0, 0, D1-D0} { Surface{s2};  }
+
+    //Compound Volume(4) = {vol2[1], vol3[1]};
+
+//Circle(12) = {11, 8, 9};
+//Line(13) = {9, 11};
+//Line Loop(14) = {12, 13};
+//Plane Surface(15) = {14};
+
+Return
+
+X = 1.;
+For r In {1:10}
+  Call Ring;
+EndFor
+
+/*
+
+Point(0) = { 0, 0, D0, lc};
+Point(1) = { R, 0, D0, lc};
+Point(2) = { 0, R, D0, lc};
+Point(3) = {-R, 0, D0, lc};
+
+// Connect up the points
+Circle(4) = {1, 0, 3};
+Line(5) = {3,1};
+Line Loop(6) = {4,5};
+Plane Surface(7) = {6};
+
+////////////////////////////////////
+// South Pole
+
+Point(8)  = { 0,  -0.0001, D0, lc};
+Point(9)  = { R,  -0.0001, D0, lc};
+Point(10) = { 0, -R, D0, lc};
+Point(11) = {-R,  -0.0001, D0, lc};
+
+// Connect up the points
+Circle(12) = {11, 8, 9};
+Line(13) = {9, 11};
+Line Loop(14) = {12, 13};
+Plane Surface(15) = {14};
+
+//////////////////////////////////////
+// Extrude magnet
+Extrude {0, 0, D1-D0} { Surface{7};  }
+Extrude {0, 0, D1-D0} { Surface{15}; }
+
+/////////////////////////////////////
+// Large Bounding box
+Point(116)  = {X0, Y0, Z0};
+Point(117)  = {X1, Y0, Z0};
+Point(118)  = {X1, Y1, Z0};
+Point(119)  = {X0, Y1, Z0};
+Line(120) = {116,117};
+Line(121) = {117,118};
+Line(122) = {118,119};
+Line(123) = {119,116};
+Line Loop(124) = {120, 121, 122, 123};
+Plane Surface(125) = {124};
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
+
+//////////////////////////////////////
+// Rings of sensitivity calculation
+//Circle(222) = {2.*1, 0, 2.*2};
+
+
+//////////////////////////////////////
+// Volumes
+
+Physical Volume("plus") = {2};
+Physical Volume("minus") = {1};
+Physical Volume("background") = {3};
+Coherence;
+*/

examples/borehole/boreholeMagnetGrid_merge.geo

@@ -0,0 +1,85 @@
+/**
+ *  Trevor P. Irons, XRI Geophysics, LLC
+ *  Test of Coloumbic Magnetic Potential
+ */
+
+lc = 1e-2;     // Target element size
+R  =  .11;     // Minimum Radius
+R2 =  .175;    // Minimum Radius
+D0 = 9.75;     // Top of mesh
+D1 = 11.25;    // Bottom of mesh
+
+////////////////////////////////////
+// North Pole
+
+Function Ring
+
+    // centre
+    p0 = newp; Point(p0) = {   0,    0, D0, lc};
+
+    // Points defining outer ring
+    p1 = newp; Point(p1) = { R,  0, D0, lc};
+    p2 = newp; Point(p2) = { 0,  R, D0, lc};
+    p3 = newp; Point(p3) = {-R,  0, D0, lc};
+    p4 = newp; Point(p4) = { 0, -R, D0, lc};
+
+    c1 = newc; Circle(c1) = {p1, p0, p2};
+    c2 = newc; Circle(c2) = {p2, p0, p3};
+    c3 = newc; Circle(c3) = {p3, p0, p4};
+    c4 = newc; Circle(c4) = {p4, p0, p1};
+
+    // Inner ring
+    p5 = newp; Point(p5) = { R2,   0, D0, lc};
+    p6 = newp; Point(p6) = {  0,  R2, D0, lc};
+    p7 = newp; Point(p7) = {-R2,   0, D0, lc};
+    p8 = newp; Point(p8) = {  0, -R2, D0, lc};
+
+    c5 = newc; Circle(c5) = {p5, p0, p6};
+    c6 = newc; Circle(c6) = {p6, p0, p7};
+    c7 = newc; Circle(c7) = {p7, p0, p8};
+    c8 = newc; Circle(c8) = {p8, p0, p5};
+
+    l1 = newl; Line Loop(l1) = {c1, c2, c3, c4, c5, c6, c7, c8};
+
+    s1 = news; Plane Surface(s1) = {l1};
+    Extrude {0,0,D1-D0} {
+      Surface{s1};
+    }
+
+Return
+
+Function Cyl
+
+    // centre
+    p0 = newp; Point(p0) = {   0,    0, D0, lc};
+
+    // Points defining outer ring
+    p1 = newp; Point(p1) = { R2,   0, D0, lc};
+    p2 = newp; Point(p2) = {  0,  R2, D0, lc};
+    p3 = newp; Point(p3) = {-R2,   0, D0, lc};
+    p4 = newp; Point(p4) = {  0, -R2, D0, lc};
+
+    c1 = newc; Circle(c1) = {p1, p0, p2};
+    c2 = newc; Circle(c2) = {p2, p0, p3};
+    c3 = newc; Circle(c3) = {p3, p0, p4};
+    c4 = newc; Circle(c4) = {p4, p0, p1};
+
+    l1 = newl; Line Loop(l1) = {c1, c2, c3, c4};
+
+    s1 = news; Plane Surface(s1) = {l1};
+    Extrude {0,0,D1-D0} {
+      Surface{s1};
+    }
+
+Return
+
+
+Call Cyl;
+/*
+For r In {1:12}
+  Call Ring;
+  R2 = R;
+  R *= 1.2;
+  lc *= 1.1;
+EndFor
+*/

examples/borehole/boreholeMagnetMinus.geo

@@ -0,0 +1,75 @@
+/**
+ *  Trevor P. Irons, XRI Geophysics, LLC
+ *  Test of Coloumbic Magnetic Potential
+ */
+
+lc = 1e-2;   // Target element size
+R = .05;   // Magnet Radius
+D0 = 10;   // Top of magnet
+D1 = 11;   // Bottom of magnet
+
+// Total Solution Space
+X0 = -2.;
+X1 =  2.;
+Y0 = -2.;
+Y1 =  2.;
+Z0 =  9.;
+Z1 = 12.;
+
+////////////////////////////////////
+// North Pole
+/*
+Point(0) = { 0, 0, D0, lc};
+Point(1) = { R, 0, D0, lc};
+Point(2) = { 0, R, D0, lc};
+Point(3) = {-R, 0, D0, lc};
+
+// Connect up the points
+Circle(4) = {1, 0, 3};
+Line(5) = {3,1};
+Line Loop(6) = {4,5};
+Plane Surface(7) = {6};
+*/
+
+
+////////////////////////////////////
+// South Pole
+
+Point(8)  = { 0,  -0.0001, D0, lc};
+Point(9)  = { R,  -0.0001, D0, lc};
+Point(10) = { 0, -R, D0, lc};
+Point(11) = {-R,  -0.0001, D0, lc};
+
+// Connect up the points
+Circle(12) = {11, 8, 9};
+Line(13) = {9, 11};
+Line Loop(14) = {12, 13};
+Plane Surface(15) = {14};
+
+
+//////////////////////////////////////
+// Extrude magnet
+//Extrude {0, 0, D1-D0} { Surface{7};  }
+Extrude {0, 0, D1-D0} { Surface{15}; }
+
+/*
+// Large Bounding box
+Point(116)  = {X0, Y0, Z0};
+Point(117)  = {X1, Y0, Z0};
+Point(118)  = {X1, Y1, Z0};
+Point(119)  = {X0, Y1, Z0};
+Line(120) = {116,117};
+Line(121) = {117,118};
+Line(122) = {118,119};
+Line(123) = {119,116};
+Line Loop(124) = {120, 121, 122, 123}; 
+Plane Surface(125) = {124};
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
+*/
+
+//////////////////////////////////////
+// Volumes
+
+Physical Volume("plus") = {2};
+Physical Volume("minus") = {1};
+//Physical Volume("background") = {3};

examples/borehole/boreholeMagnetPlus.geo

@@ -0,0 +1,175 @@
+/**
+ *  Trevor P. Irons, XRI Geophysics, LLC
+ *  Test of Coloumbic Magnetic Potential
+ */
+
+
+Function Cyl
+    // centre
+    pp0 = newp; Point(pp0) = {   0,    0, D0, lc2};
+
+    // Points defining outer ring
+    pp1 = newp; Point(pp1) = { R2,   0, D0, lc2};
+    pp2 = newp; Point(pp2) = {  0,  R2, D0, lc2};
+    pp3 = newp; Point(pp3) = {-R2,   0, D0, lc2};
+    pp4 = newp; Point(pp4) = {  0, -R2, D0, lc2};
+
+    cc1 = newc; Circle(cc1) = {pp1, pp0, pp2};
+    cc2 = newc; Circle(cc2) = {pp2, pp0, pp3};
+    cc3 = newc; Circle(cc3) = {pp3, pp0, pp4};
+    cc4 = newc; Circle(cc4) = {pp4, pp0, pp1};
+
+    ll1 = newl; Line Loop(ll1) = {cc1, cc2, cc3, cc4};
+
+    ss1 = news; Plane Surface(ss1) = {ll1};
+    //v1 = newv; v1 =
+    //newv vol2;
+    vol2[] = Extrude {0, 0, D1-D0} { Surface{ss1}; };
+
+Return
+
+lc = 5e-2;   // Target element size
+R = .05;   // Magnet Radius
+D0 = 10;   // Top of magnet
+D1 = 11;   // Bottom of magnet
+
+
+p0 = newp; Point(p0) = { 0, 0, D0, lc};
+p1 = newp; Point(p1) = { R, 0, D0, lc};
+p2 = newp; Point(p2) = { 0, R, D0, lc};
+p3 = newp; Point(p3) = {-R, 0, D0, lc};
+p4 = newp; Point(p4) = { 0,-R, D0, lc};
+
+c1 = newc; Circle(c1) = {p1, p0, p2};
+c2 = newc; Circle(c2) = {p2, p0, p3};
+c3 = newc; Circle(c3) = {p3, p0, p4};
+c4 = newc; Circle(c4) = {p4, p0, p1};
+
+
+ll0 = newl; Line Loop(ll0) = {1,2,3,4};
+Plane Surface(6) = {5};
+vol2[] = Extrude {0,0,D1-D0} {
+  Surface{6};
+};
+
+
+/*
+vol1[] = Extrude {0, 0, D1-D0} { Line{c1:c4}; };
+*/
+////////////////////////////////////////////////////
+// Big Cylinder
+////////////////////////////////////////////////////
+/*
+lc2 = 1e-2;    // Target element size
+R2 =  .175;    // Minimum Radius
+D0 = 9.75;     // Top of mesh
+D1 = 11.25;    // Bottom of mesh
+
+Call Cyl;
+
+//newv vol3;
+//Compound Volume(vol3) = {vol1[1], vol2[1]};
+
+
+Field[1] = Attractor;
+Field[1].NodesList = {p0,p1,p2,p3,p4};
+
+Field[2] = Threshold;
+Field[2].IField = 1;
+Field[2].LcMin = lc / 5;
+Field[2].LcMax = lc;
+Field[2].DistMin = 0.05;
+Field[2].DistMax = 0.5;
+
+// Use minimum of all the fields as the background field
+Field[3] = Min;
+Field[3].FieldsList = {2};
+Background Field = 2;
+
+// Don't extend the elements sizes from the boundary inside the domain
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromPoints = 1 ;
+*/
+Mesh 3;
+
+/*
+Delete { Volume{vol1[1]}; }
+Delete { Volume{vol2[1]}; }
+
+Delete { Point{p0}; }
+Delete { Point{p1}; }
+Delete { Point{p2}; }
+Delete { Point{p3}; }
+Delete { Point{p4}; }
+
+Delete { Line{c1}; }
+Delete { Line{c2}; }
+Delete { Line{c3}; }
+Delete { Line{c4}; }
+*/
+//Delete { Line{ll0}; }
+//Delete { Volume{vol2}; }
+
+
+
+
+
+
+/*
+////////////////////////////////////
+// South Pole
+
+Point(8)  = { 0,  -0.0001, D0, lc};
+Point(9)  = { R,  -0.0001, D0, lc};
+Point(10) = { 0, -R, D0, lc};
+Point(11) = {-R,  -0.0001, D0, lc};
+
+// Connect up the points
+Circle(12) = {11, 8, 9};
+Line(13) = {9, 11};
+Line Loop(14) = {12, 13};
+Plane Surface(15) = {14};
+*/
+
+//////////////////////////////////////
+// Extrude magnet
+//Extrude {0, 0, D1-D0} { Surface{7};  }
+//Extrude {0, 0, D1-D0} { Surface{15}; }
+
+/*
+// Large Bounding box
+Point(116)  = {X0, Y0, Z0};
+Point(117)  = {X1, Y0, Z0};
+Point(118)  = {X1, Y1, Z0};
+Point(119)  = {X0, Y1, Z0};
+Line(120) = {116,117};
+Line(121) = {117,118};
+Line(122) = {118,119};
+Line(123) = {119,116};
+Line Loop(124) = {120, 121, 122, 123};
+Plane Surface(125) = {124};
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
+*/
+
+//////////////////////////////////////
+// Volumes
+
+//Physical Volume("plus") = {2};
+//Physical Volume("minus") = {1};
+//Physical Volume("background") = {3};
+/*
+Delete {
+  Surface{9};
+}
+Delete {
+  Surface{13};
+}
+Delete {
+  Surface{17};
+}
+Delete {
+  Surface{21};
+}
+*/
+Compound Volume(49) = {1};
+Compound Volume(50) = {49, 1};

examples/borehole/sphere.geo

@@ -0,0 +1,146 @@
+/* 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      08/08/2014 12:19:20 PM
+ * @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
+ */
+
+
+D0 = 10;          // Top of magnet
+D1 = 11;          // Bottom of magnet
+radius = 2.25;     // Radius of the damn thing
+
+lc = radius*10;   //  0.25;   // Target element size
+
+// Total Solution Space
+Box = 10*radius; // The down side of potential
+X0 = -Box;
+X1 =  Box;
+Y0 = -Box;
+Y1 =  Box;
+Z0 = -Box;
+Z1 =  Box;
+
+cellSize=lc; //300;
+dd = 0 ; //  1e-5; //cellSize; // .01;
+pio2=Pi/2;
+
+// Calculate offset effect
+theta = Asin(dd/radius);
+rr = radius * Cos(theta);
+
+///////////////////////////////////
+// Positive half sphere
+// create inner 1/8 shell
+p0 = newp;
+Point(p0)   = {      0,   0,       0, cellSize}; // origin
+Point(p0+1) = {    -rr,   0,      dd, cellSize};
+Point(p0+2) = {      0,  rr,      dd, cellSize};
+Point(p0+3) = {      0,   0,  radius, cellSize};
+Point(p0+4) = {      0,   0,      dd, cellSize}; // origin
+
+c0 = newc;
+Circle(c0  ) = {p0+1, p0+4, p0+2};       // Tricky, This one needs to be offset!
+Circle(c0+1) = {p0+3, p0, p0+1};
+Circle(c0+2) = {p0+3, p0, p0+2};
+
+Line Loop(10) = {c0, -(c0+2), c0+1} ;
+Ruled Surface (60) = {10};
+
+////////////////////////////////////////////////////////////
+// Negative half sphere
+p = newp;
+Point(  p) = {      0,      0,            0, cellSize};
+Point(p+1) = {    -rr,      0,          -dd, cellSize};
+Point(p+2) = {      0,     rr,          -dd, cellSize};
+Point(p+3) = {      0,      0,      -radius, cellSize};
+Point(p+4) = {      0,      0,          -dd, cellSize};
+
+cc = newc;
+Circle(cc  ) = {p+1, p+4, p+2};
+Circle(cc+1) = {p+3, p,   p+1};
+Circle(cc+2) = {p+3, p,   p+2};
+
+Circle(cc+3) = {p+3, p, p+2};
+Circle(cc+4) = {p+3, p, p+2};
+Circle(cc+5) = {p+3, p, p+2};
+
+ccl = newl;
+Line(ccl) = { p0+3, p+3 };
+
+Line Loop(11) = {cc, -(cc+2), cc+1} ;
+Ruled Surface (61) = {11};
+
+// create remaining 7/8 inner shells
+t1[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{60};}};
+t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
+t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
+//
+t4[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{61};}};
+t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
+t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
+
+/////////////////////////////////////
+// Large Bounding box
+pp = newp;
+Point(pp)    = {X0, Y0, Z0, lc};
+Point(pp+1)  = {X1, Y0, Z0, lc};
+Point(pp+2)  = {X1, Y1, Z0, lc};
+Point(pp+3)  = {X0, Y1, Z0, lc};
+
+
+lv = newl;
+Line(lv) = {pp,pp+1};
+Line(lv+1) = {pp+1,pp+2};
+Line(lv+2) = {pp+2,pp+3};
+Line(lv+3) = {pp+3,pp};
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
+
+// Hard coded doom
+Plane Surface(125) = {lv+4};
+
+//v = newv;
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
+
+/* This is GOOD */
+Surface{60} In Volume{v[1]};
+Surface{t1[0]} In Volume{v[1]};
+Surface{t2[0]} In Volume{v[1]};
+Surface{t3[0]} In Volume{v[1]};
+
+Surface{61} In Volume{v[1]};
+Surface{t4[0]} In Volume{v[1]};
+Surface{t5[0]} In Volume{v[1]};
+Surface{t6[0]} In Volume{v[1]};
+
+///////////////////////////////////////////////
+// Attractor Field
+
+Field[1] = Attractor;
+Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
+
+Field[2] = MathEval;
+//Field[2].F = Sprintf("(.25 - F1)^2 + %g", cellSize );  // WORKS
+Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, cellSize/2. );
+Background Field = 2;
+
+// Don't extend the elements sizes from the boundary inside the domain
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
+
+Physical Volume(1) = {v[1]};
+
+// To create the mesh run
+// gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
+//gmsh -3 -format msh1 -o outfile.msh sphere.geo

examples/borehole/sphere.sh

@@ -0,0 +1,26 @@
+
+#gmsh -3 -format vtk -o sphere.vtk  sphere.geo
+#../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
+#../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 
+
+# Gravity problem
+../utVTKEdgeGsphere sphere.vtk .25 gravity  sphereGrav.vtu
+../utFEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
+#paraview --state=sliceSphere.pvsm
+
+
+# --data=sphereOut.vtu
+
+# even finer meshing
+# radius = 4, u = \pm .9              analytic = \pm 1.33
+
+# Finer meshing
+# radius = 5, u = \pm 1.65            analytic = \pm 1.67
+# radius = 4, u = \pm 1.09            analytic = \pm 1.33
+
+# radius = 5, u = \pm 2.5             analytic = \pm 1.67
+# radius = 4, u = \pm 1.66            analytic = \pm 1.33
+# radius = 3, u= \pm .926 -.985       analytic = \pm 1
+# radius = 2, u = \pm.42              analytic = \pm .667
+# for radius=1,  u = \pm .111         analytic = \pm .333 
+# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167

examples/borehole/untitled.geo

@@ -0,0 +1 @@
+Point(1) = {-0.25, 0, 10, 1.0};