Selaa lähdekoodia

Merge branch 'iss0'

iss2
Trevor Irons 8 vuotta sitten
vanhempi
commit
2b3ccbc334

+ 6
- 2
CMakeLists.txt Näytä tiedosto

@@ -6,6 +6,10 @@ add_subdirectory("src")
6 6
 
7 7
 add_library( fem4ellipticpde ${FEMSOURCE} )  
8 8
 target_link_libraries(fem4ellipticpde "lemmacore")
9
-else()
10
-	message( STATUS "FEM4EllipticPDE requires Lemma to be build with VTK" )
9
+
10
+install ( TARGETS fem4ellipticpde DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
11
+
12
+if (LEMMABUILDEXAMPLES)
13
+	add_subdirectory(examples)
14
+endif()
11 15
 endif()

+ 33
- 0
examples/CMakeLists.txt Näytä tiedosto

@@ -0,0 +1,33 @@
1
+add_executable( FEM4EllipticPDE_bhmag FEM4EllipticPDE_bhmag.cpp  )
2
+target_link_libraries(  FEM4EllipticPDE_bhmag  "lemmacore" "fem4ellipticpde")
3
+
4
+add_executable( FEM4EllipticPDE FEM4EllipticPDE.cpp  )
5
+target_link_libraries(  FEM4EllipticPDE  "lemmacore" "fem4ellipticpde")
6
+
7
+add_executable( merge merge.cpp  )
8
+target_link_libraries(  merge  "lemmacore" "fem4ellipticpde")
9
+
10
+add_executable( VTKDC VTKDC.cpp  )
11
+target_link_libraries( VTKDC  "lemmacore" "fem4ellipticpde")
12
+
13
+add_executable( VTKEdgeG VTKEdgeG.cpp  )
14
+target_link_libraries( VTKEdgeG  "lemmacore" "fem4ellipticpde")
15
+
16
+add_executable( VTKEdgeGsphere VTKEdgeGsphere.cpp  )
17
+target_link_libraries( VTKEdgeGsphere  "lemmacore" "fem4ellipticpde")
18
+
19
+#INSTALL_TARGETS( "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
20
+INSTALL_TARGETS( "/share/FEM4EllipticPDE/"
21
+	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKEdgeGsphere
22
+)
23
+
24
+install( DIRECTORY "borehole"
25
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
26
+	PATTERN .svn EXCLUDE
27
+)
28
+
29
+install( DIRECTORY "DC"
30
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
31
+	PATTERN .svn EXCLUDE
32
+)
33
+

+ 126
- 0
examples/DC/DC.yaml Näytä tiedosto

@@ -0,0 +1,126 @@
1
+!<DCSurvey>
2
+NumberOfElectrodes: 9
3
+Injections:
4
+  J-0:
5
+    A: !<DCIPElectrode> &1
6
+      Node_ID: -1
7
+      Label: L0-E0
8
+      Location: !<Vector3r>
9
+        -
10
+          - 13.23
11
+          - 18.1
12
+          - -25
13
+    B: !<DCIPElectrode> &2
14
+      Location: !<Vector3r>
15
+        -
16
+          - 13.23
17
+          - 18.1
18
+          - -25
19
+      Node_ID: -1
20
+      Label: L0-E1
21
+    J: 1
22
+    Measurements:
23
+      Number: 4
24
+      V-0:
25
+        M: !<DCIPElectrode> &3
26
+          Location: !<Vector3r>
27
+            -
28
+              - 13.23
29
+              - 18.1
30
+              - -25
31
+          Node_ID: -1
32
+          Label: L0-E2
33
+        N: !<DCIPElectrode> &4
34
+          Location: !<Vector3r>
35
+            -
36
+              - 13.23
37
+              - 18.1
38
+              - -25
39
+          Node_ID: -1
40
+          Label: L0-E3
41
+      V-1:
42
+        M: !<DCIPElectrode> &5
43
+          Location: !<Vector3r>
44
+            -
45
+              - 13.23
46
+              - 18.1
47
+              - -25
48
+          Node_ID: -1
49
+          Label: L0-E4
50
+        N: !<DCIPElectrode> &6
51
+          Location: !<Vector3r>
52
+            -
53
+              - 13.23
54
+              - 18.1
55
+              - -25
56
+          Node_ID: -1
57
+          Label: L0-E5
58
+      V-2:
59
+        M: *6
60
+        N: !<DCIPElectrode> &7
61
+          Location: !<Vector3r>
62
+            -
63
+              - 13.23
64
+              - 18.1
65
+              - -25
66
+          Node_ID: -1
67
+          Label: L0-E6
68
+      V-3:
69
+        M: !<DCIPElectrode> &8
70
+          Location: !<Vector3r>
71
+            -
72
+              - 13.23
73
+              - 18.1
74
+              - -25
75
+          Node_ID: -1
76
+          Label: L0-E7
77
+        N: !<DCIPElectrode> &9
78
+          Location: !<Vector3r>
79
+            -
80
+              - 13.23
81
+              - 18.1
82
+              - -25
83
+          Node_ID: -1
84
+          Label: L0-E8
85
+  J-1:
86
+    A: *2
87
+    B: *3
88
+    J: 1
89
+    Measurements:
90
+      Number: 3
91
+      V-0:
92
+        M: *3
93
+        N: *4
94
+      V-1:
95
+        M: *5
96
+        N: *6
97
+      V-2:
98
+        M: *6
99
+        N: *7
100
+  J-2:
101
+    A: *3
102
+    B: *4
103
+    J: 1
104
+    Measurements:
105
+      Number: 3
106
+      V-0:
107
+        M: *3
108
+        N: *4
109
+      V-1:
110
+        M: *5
111
+        N: *6
112
+      V-2:
113
+        M: *6
114
+        N: *7
115
+Electrodes:
116
+  L0-E2: *3
117
+  L0-E3: *4
118
+  L0-E4: *5
119
+  L0-E5: *6
120
+  L0-E6: *7
121
+  L0-E7: *8
122
+  L0-E8: *9
123
+  L0-E0: *1
124
+  L0-E1: *2
125
+---
126
+

+ 198
- 0
examples/DC/tunnel.geo Näytä tiedosto

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

+ 493
- 0
examples/FEM4EllipticPDE.cpp Näytä tiedosto

@@ -0,0 +1,493 @@
1
+// ===========================================================================
2
+//
3
+//       Filename:  utFEM4EllipticPDE.cpp
4
+//
5
+//        Created:  08/16/12 19:49:10
6
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
7
+//
8
+//         Author:  Trevor Irons (ti)
9
+//
10
+//   Organisation:  Colorado School of Mines (CSM)
11
+//                  United States Geological Survey (USGS)
12
+//
13
+//          Email:  tirons@mines.edu, tirons@usgs.gov
14
+//
15
+//  This program is free software: you can redistribute it and/or modify
16
+//  it under the terms of the GNU General Public License as published by
17
+//  the Free Software Foundation, either version 3 of the License, or
18
+//  (at your option) any later version.
19
+//
20
+//  This program is distributed in the hope that it will be useful,
21
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
22
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23
+//  GNU General Public License for more details.
24
+//
25
+//  You should have received a copy of the GNU General Public License
26
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
27
+//
28
+// ===========================================================================
29
+
30
+/**
31
+  @file
32
+  @author   Trevor Irons
33
+  @date     08/16/12
34
+  @version   0.0
35
+ **/
36
+
37
+#include "FEM4EllipticPDE.h"
38
+#include "vtkSphere.h"
39
+#include "vtkRectilinearGrid.h"
40
+#include "vtkFloatArray.h"
41
+#include <vtkIdList.h>
42
+#include <vtkIdTypeArray.h>
43
+#include <vtkCell.h>
44
+#include <vtkTriangleFilter.h>
45
+#include <vtkDataSetSurfaceFilter.h>
46
+#include <vtkExtractCells.h>
47
+#include <vtkGeometryFilter.h>
48
+#include <vtkUnstructuredGrid.h>
49
+#include <vtkExtractEdges.h>
50
+#include <vtkImplicitDataSet.h>
51
+
52
+// debugging stuff
53
+#include <vtkDataSetMapper.h>
54
+#include <vtkSelectionNode.h>
55
+#include <vtkExtractSelection.h>
56
+#include <vtkSelection.h>
57
+#include <vtkVertexGlyphFilter.h>
58
+
59
+#define RENDERTEST
60
+//#define CONNECTTEST
61
+
62
+#ifdef RENDERTEST
63
+#include "vtkRectilinearGridGeometryFilter.h"
64
+#include "vtkRectilinearGridOutlineFilter.h"
65
+#include "vtkPolyDataMapper.h"
66
+#include "vtkActor.h"
67
+#include "vtkRenderWindowInteractor.h"
68
+#include "vtkRenderer.h"
69
+#include "vtkRenderWindow.h"
70
+#include "vtkCamera.h"
71
+#include "vtkProperty.h"
72
+#include <vtkDataSetMapper.h>
73
+#include <vtkSelectionNode.h>
74
+#include <vtkSelection.h>
75
+#include <vtkExtractSelection.h>
76
+#include <vtkExtractEdges.h>
77
+#include <vtkVertexGlyphFilter.h>
78
+#include <vtkTriangleFilter.h>
79
+#include <vtkImplicitHalo.h>
80
+#endif
81
+
82
+#ifdef CONNECTTEST
83
+#include <vtkVersion.h>
84
+#include <vtkSmartPointer.h>
85
+#include <vtkPolyData.h>
86
+#include <vtkCellData.h>
87
+#include <vtkDoubleArray.h>
88
+#include <vtkDataSet.h>
89
+#include <vtkSphereSource.h>
90
+#include <vtkTriangleFilter.h>
91
+#include <vtkExtractEdges.h>
92
+#include <vtkDataSetMapper.h>
93
+#include <vtkActor.h>
94
+#include <vtkRenderWindow.h>
95
+#include <vtkRenderer.h>
96
+#include <vtkRenderWindowInteractor.h>
97
+#include <vtkSelectionNode.h>
98
+#include <vtkSelection.h>
99
+#include <vtkExtractSelection.h>
100
+#include <vtkProperty.h>
101
+#include <vtkVertexGlyphFilter.h>
102
+#endif
103
+
104
+vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id);
105
+vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id);
106
+
107
+// This function is copyright (C) 2012 Joeseph Cappriotti
108
+vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid);
109
+
110
+using namespace Lemma;
111
+
112
+int main() {
113
+
114
+    // Define Sigma term
115
+    // Defaults to 1, (Poisson equation)
116
+
117
+    // Define source (G) term
118
+    //vtkSphere *Sphere = vtkSphere::New();
119
+    vtkImplicitHalo *Sphere = vtkImplicitHalo::New();
120
+        Sphere->SetCenter (0, 0, 0);
121
+        Sphere->SetRadius (60);
122
+        Sphere->SetFadeOut(.95);
123
+
124
+    ////////////////////////////////////////////
125
+    // Define solution mesh
126
+
127
+    // Create a rectilinear grid by defining three arrays specifying the
128
+    // coordinates in the x-y-z directions.
129
+
130
+    int    nx = 160;
131
+    double dx = 2.5;
132
+    double ox = -200;
133
+    vtkFloatArray *xCoords = vtkFloatArray::New();
134
+    for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
135
+
136
+    int    ny = 160;
137
+    double dy = 2.5;
138
+    double oy = -200;
139
+    vtkFloatArray *yCoords = vtkFloatArray::New();
140
+    for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
141
+
142
+    int    nz = 160;
143
+    double dz = 2.5;
144
+    double oz = -200;
145
+    vtkFloatArray *zCoords = vtkFloatArray::New();
146
+    for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
147
+
148
+    vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
149
+        rGrid->SetDimensions(nx+1, ny+1, nz+1);
150
+        rGrid->SetExtent(0, nx, 0, ny, 0,  nz);
151
+        rGrid->SetXCoordinates(xCoords);
152
+        rGrid->SetYCoordinates(yCoords);
153
+        rGrid->SetZCoordinates(zCoords);
154
+
155
+    vtkRectilinearGrid *sigGrid = vtkRectilinearGrid::New();
156
+        sigGrid->SetDimensions(nx+1, ny+1, nz+1);
157
+        sigGrid->SetExtent(0, nx, 0, ny, 0, nz);
158
+        sigGrid->SetXCoordinates(xCoords);
159
+        sigGrid->SetYCoordinates(yCoords);
160
+        sigGrid->SetZCoordinates(zCoords);
161
+        vtkDoubleArray *sigArray = vtkDoubleArray::New();
162
+        sigArray->SetNumberOfComponents(1);
163
+        for (int ii=0; ii<(nx+1)*(ny+1)*(nz+1); ++ii) {
164
+            sigArray->InsertTuple1(ii, 1.);
165
+        }
166
+        sigArray->SetName("sigma");
167
+        sigGrid->GetPointData()->AddArray(sigArray);
168
+        sigGrid->GetPointData()->SetScalars(sigArray);
169
+
170
+    vtkImplicitDataSet* implSigma = vtkImplicitDataSet::New();
171
+        implSigma->SetDataSet(sigGrid);
172
+        implSigma->SetOutValue(0.);
173
+        implSigma->SetOutGradient(0., 0., 0.);
174
+
175
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
176
+        Solver->SetGFunction(Sphere);
177
+        //Solver->SetSigmaFunction(implSigma);
178
+        Solver->SetBoundaryStep(.5);
179
+        Solver->SetGrid(rGrid);
180
+        Solver->Solve("FEM_results.vtr");
181
+
182
+    vtkXMLDataSetWriter *Writer = vtkXMLDataSetWriter::New();
183
+        Writer->SetInputData(sigGrid);
184
+        std::string fname("sigma.vtr");
185
+        //Writer->Update();
186
+        //fname.append( Writer->GetDefaultFileExtension() );
187
+        Writer->SetFileName(fname.c_str());
188
+        Writer->Write();
189
+        Writer->Delete();
190
+
191
+    // Clean up
192
+    Sphere->Delete();
193
+    Solver->Delete();
194
+    rGrid->Delete();
195
+    xCoords->Delete();
196
+    yCoords->Delete();
197
+    zCoords->Delete();
198
+
199
+}
200
+
201
+// This function is copyright (C) 2012 Joeseph Cappriotti
202
+vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid){
203
+    vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
204
+    vtkSmartPointer<vtkIdList> cellList = vtkSmartPointer<vtkIdList>::New();
205
+    grid->GetPointCells(id0, cellList);
206
+    for(int i=0;i<cellList->GetNumberOfIds();++i){
207
+        vtkCell* cell = grid->GetCell(cellList->GetId(i));
208
+        if(cell->GetNumberOfEdges() > 0){
209
+            for(int j=0; j<cell->GetNumberOfEdges(); ++j){
210
+                vtkCell* edge = cell->GetEdge(j);
211
+                vtkIdList* edgePoints=edge->GetPointIds();
212
+                if(edgePoints->GetId(0)==id0){
213
+                    pointIds->InsertUniqueId(edgePoints->GetId(1));
214
+                }else if(edgePoints->GetId(1)==id0){
215
+                    pointIds->InsertUniqueId(edgePoints->GetId(0));
216
+                }
217
+            }
218
+        }
219
+    }
220
+    return pointIds;
221
+}
222
+
223
+vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id) {
224
+
225
+    std::cout << "number of points " << mesh->GetNumberOfPoints() << std::endl;
226
+    std::cout << "number of cells " << mesh->GetNumberOfCells() << std::endl;
227
+
228
+    vtkSmartPointer<vtkIdList> connectedVertices = vtkSmartPointer<vtkIdList>::New();
229
+
230
+    //get all cells that vertex 'id' is a part of
231
+    vtkSmartPointer<vtkIdList> cellIdList = vtkSmartPointer<vtkIdList>::New();
232
+    mesh->GetPointCells(id, cellIdList);   // cells using point
233
+    //mesh->GetCellPoints(id, cellIdList); // points defining cell
234
+
235
+    cout << "Vertex " << id << " is used in cell(s) ";
236
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
237
+        cout << cellIdList->GetId(i) << ", ";
238
+    }
239
+    cout << endl;
240
+
241
+    // TODO this is where things don't work
242
+    // for each cell, figure out what vertices are connected
243
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
244
+
245
+        cout << "\tcell id " << i << " : " << cellIdList->GetId(i) << endl;
246
+        std::cout << "\tcell has " << mesh->GetCell(i)->GetNumberOfPoints() << " points\n";
247
+
248
+        vtkExtractCells *Cell = vtkExtractCells::New();
249
+            Cell->SetInputData(mesh);
250
+            Cell->AddCellRange(i, i);
251
+            Cell->Update();
252
+
253
+        //std::cout << *Cell->GetOutput() << std::endl;
254
+
255
+        // extract surface
256
+        vtkDataSetSurfaceFilter *surf = vtkDataSetSurfaceFilter::New();
257
+            surf->UseStripsOn();
258
+            surf->SetInputData(Cell->GetOutput());
259
+            surf->Update();
260
+
261
+        //std::cout << *surf->GetOutput() << std::endl;
262
+
263
+        vtkSmartPointer<vtkTriangleFilter> triangleFilter =
264
+             vtkSmartPointer<vtkTriangleFilter>::New();
265
+        triangleFilter->SetInputData( surf->GetOutput() );
266
+            triangleFilter->Update();
267
+
268
+        vtkSmartPointer<vtkExtractEdges> extractEdges =
269
+            vtkSmartPointer<vtkExtractEdges>::New();
270
+        extractEdges->SetInputConnection(triangleFilter->GetOutputPort());
271
+            extractEdges->Update();
272
+
273
+        vtkSmartPointer<vtkPolyData> pmesh = extractEdges->GetOutput();
274
+        //vtkSmartPointer<vtkPolyData> pmesh = triangleFilter->GetOutput();
275
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
276
+//
277
+//         //vtkSmartPointer<vtkIdTypeArray> ids =
278
+//         //vtkSmartPointer<vtkIdTypeArray>::New();
279
+//         //ids->SetNumberOfComponents(1);
280
+//
281
+        vtkSmartPointer<vtkIdList> pointIdList = vtkSmartPointer<vtkIdList>::New();
282
+        //mesh->GetPointCells(cellIdList->GetId(i), pointIdList); // returns cells using point
283
+        mesh->GetCellPoints(cellIdList->GetId(i), pointIdList);   // points defining cell
284
+
285
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(pmesh, id);
286
+        vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
287
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(triangleFilter->GetOutput(), 1);
288
+        std::cout << "\tPoly Connected vertices: ";
289
+        for(vtkIdType ip = 0; ip < pconnectedVertices->GetNumberOfIds(); ip++) {
290
+            std::cout << pconnectedVertices->GetId(ip) << " ";
291
+            //ids->InsertNextValue(connectedVertices->GetId(i));
292
+            connectedVertices->InsertNextId( pointIdList->GetId( pconnectedVertices->GetId(ip) ) );
293
+        }
294
+        std::cout << std::endl;
295
+
296
+
297
+        cout << "\tcell " << i << " contains these points ";
298
+        for(vtkIdType ii = 0; ii < pointIdList->GetNumberOfIds(); ii++) {
299
+            cout << pointIdList->GetId(ii) << ", ";
300
+        }
301
+        cout << endl;
302
+
303
+        /*
304
+        cout << "\tEnd points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
305
+
306
+        if(pointIdList->GetId(0) != id) {
307
+            //cout << "Connected to " << pointIdList->GetId(0) << endl;
308
+            connectedVertices->InsertNextId(pointIdList->GetId(0));
309
+        } else {
310
+            //cout << "Connected to " << pointIdList->GetId(1) << endl;
311
+            connectedVertices->InsertNextId(pointIdList->GetId(1));
312
+        }
313
+        */
314
+
315
+
316
+        #ifdef RENDERTEST
317
+        vtkSmartPointer<vtkDataSetMapper> sphereMapper =
318
+        vtkSmartPointer<vtkDataSetMapper>::New();
319
+            sphereMapper->SetInputConnection(extractEdges->GetOutputPort());
320
+        vtkSmartPointer<vtkActor> sphereActor =
321
+            vtkSmartPointer<vtkActor>::New();
322
+            sphereActor->SetMapper(sphereMapper);
323
+
324
+
325
+        vtkSmartPointer<vtkIdTypeArray> ids =
326
+            vtkSmartPointer<vtkIdTypeArray>::New();
327
+        ids->SetNumberOfComponents(1);
328
+
329
+        std::cout << "RENDERTEST Connected vertices: ";
330
+        for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); i++) {
331
+            std::cout << connectedVertices->GetId(i) << " ";
332
+            ids->InsertNextValue(connectedVertices->GetId(i));
333
+        }
334
+        std::cout << std::endl;
335
+
336
+        vtkSmartPointer<vtkDataSetMapper> connectedVertexMapper =
337
+        vtkSmartPointer<vtkDataSetMapper>::New();
338
+
339
+        {
340
+            vtkSmartPointer<vtkSelectionNode> selectionNode =
341
+                vtkSmartPointer<vtkSelectionNode>::New();
342
+            selectionNode->SetFieldType(vtkSelectionNode::POINT);
343
+            selectionNode->SetContentType(vtkSelectionNode::INDICES);
344
+            selectionNode->SetSelectionList(ids);
345
+
346
+            vtkSmartPointer<vtkSelection> selection =
347
+                vtkSmartPointer<vtkSelection>::New();
348
+            selection->AddNode(selectionNode);
349
+
350
+            vtkSmartPointer<vtkExtractSelection> extractSelection =
351
+                vtkSmartPointer<vtkExtractSelection>::New();
352
+
353
+            extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
354
+            extractSelection->SetInputData(1, selection);
355
+            extractSelection->Update();
356
+
357
+            vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
358
+                vtkSmartPointer<vtkVertexGlyphFilter>::New();
359
+            glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
360
+            glyphFilter->Update();
361
+
362
+            connectedVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
363
+        }
364
+
365
+        vtkSmartPointer<vtkActor> connectedVertexActor =
366
+            vtkSmartPointer<vtkActor>::New();
367
+        connectedVertexActor->SetMapper(connectedVertexMapper);
368
+        connectedVertexActor->GetProperty()->SetColor(1,0,0);
369
+        connectedVertexActor->GetProperty()->SetPointSize(5);
370
+
371
+        vtkSmartPointer<vtkDataSetMapper> queryVertexMapper =
372
+            vtkSmartPointer<vtkDataSetMapper>::New();
373
+
374
+        {
375
+            vtkSmartPointer<vtkIdTypeArray> ids =
376
+                vtkSmartPointer<vtkIdTypeArray>::New();
377
+            ids->SetNumberOfComponents(1);
378
+            ids->InsertNextValue(id);
379
+
380
+            vtkSmartPointer<vtkSelectionNode> selectionNode =
381
+                vtkSmartPointer<vtkSelectionNode>::New();
382
+            selectionNode->SetFieldType(vtkSelectionNode::POINT);
383
+            selectionNode->SetContentType(vtkSelectionNode::INDICES);
384
+            selectionNode->SetSelectionList(ids);
385
+
386
+            vtkSmartPointer<vtkSelection> selection =
387
+                vtkSmartPointer<vtkSelection>::New();
388
+            selection->AddNode(selectionNode);
389
+
390
+            vtkSmartPointer<vtkExtractSelection> extractSelection =
391
+                vtkSmartPointer<vtkExtractSelection>::New();
392
+
393
+            extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
394
+            extractSelection->SetInputData(1, selection);
395
+            extractSelection->Update();
396
+
397
+            vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
398
+                vtkSmartPointer<vtkVertexGlyphFilter>::New();
399
+            glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
400
+            glyphFilter->Update();
401
+
402
+            queryVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
403
+        }
404
+
405
+        vtkSmartPointer<vtkActor> queryVertexActor =
406
+            vtkSmartPointer<vtkActor>::New();
407
+        queryVertexActor->SetMapper(queryVertexMapper);
408
+        queryVertexActor->GetProperty()->SetColor(0,1,0);
409
+        queryVertexActor->GetProperty()->SetPointSize(5);
410
+
411
+        // Create the usual rendering stuff.
412
+        vtkRenderer *renderer = vtkRenderer::New();
413
+        vtkRenderWindow *renWin = vtkRenderWindow::New();
414
+        renWin->AddRenderer(renderer);
415
+        vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
416
+        iren->SetRenderWindow(renWin);
417
+
418
+        renderer->AddActor(queryVertexActor);
419
+        renderer->AddActor(connectedVertexActor);
420
+        renderer->AddActor(sphereActor);
421
+
422
+        //renderer->AddActor(wireActor);
423
+        //renderer->AddActor(owireActor);
424
+        renderer->SetBackground(.3,.2,.1);
425
+        renderer->ResetCamera();
426
+        renderer->GetActiveCamera()->Elevation(-150.0);
427
+        renderer->GetActiveCamera()->Azimuth(0.0);
428
+        renderer->GetActiveCamera()->Roll(0.0);
429
+        //renderer->GetActiveCamera()->Pitch(-45.0);
430
+        renderer->GetActiveCamera()->Zoom(1.0);
431
+
432
+        renWin->SetSize(300,300);
433
+
434
+        // interact with data
435
+        renWin->Render();
436
+        iren->Start();
437
+
438
+        // clean up rendering stuff
439
+        //plane->Delete();
440
+        //rgridMapper->Delete();
441
+        //wireActor->Delete();
442
+        //owireActor->Delete();
443
+        //outline->Delete();
444
+        //outlineMapper->Delete();
445
+        renderer->Delete();
446
+        renWin->Delete();
447
+        iren->Delete();
448
+        #endif // RENDERTEST
449
+
450
+        surf->Delete();
451
+        Cell->Delete();
452
+    }
453
+
454
+    return connectedVertices;
455
+}
456
+
457
+// from http://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/VertexConnectivity
458
+vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id) {
459
+
460
+    std::cout << "mesh points " << mesh->GetNumberOfPoints() << std::endl;
461
+
462
+    vtkSmartPointer<vtkIdList> connectedVertices =
463
+        vtkSmartPointer<vtkIdList>::New();
464
+
465
+    //get all cells that vertex 'id' is a part of
466
+    vtkSmartPointer<vtkIdList> cellIdList =
467
+        vtkSmartPointer<vtkIdList>::New();
468
+    mesh->GetPointCells(id, cellIdList);
469
+
470
+    /*
471
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
472
+        cout << cellIdList->GetId(i) << ", ";
473
+    }
474
+    cout << endl;
475
+    */
476
+
477
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
478
+        //cout << "id " << i << " : " << cellIdList->GetId(i) << endl;
479
+
480
+        vtkSmartPointer<vtkIdList> pointIdList =
481
+            vtkSmartPointer<vtkIdList>::New();
482
+        mesh->GetCellPoints(cellIdList->GetId(i), pointIdList);
483
+        //cout << "End points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
484
+        if(pointIdList->GetId(0) != id) {
485
+            //cout << "Connected to " << pointIdList->GetId(0) << endl;
486
+            connectedVertices->InsertNextId(pointIdList->GetId(0));
487
+        } else {
488
+            //cout << "Connected to " << pointIdList->GetId(1) << endl;
489
+            connectedVertices->InsertNextId(pointIdList->GetId(1));
490
+        }
491
+    }
492
+    return connectedVertices;
493
+}

+ 173
- 0
examples/FEM4EllipticPDE_bhmag.cpp Näytä tiedosto

@@ -0,0 +1,173 @@
1
+// ===========================================================================
2
+//
3
+//       Filename:  utFEM4EllipticPDE.cpp
4
+//
5
+//        Created:  08/16/12 19:49:10
6
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
7
+//
8
+//         Author:  Trevor Irons (ti)
9
+//
10
+//   Organisation:  Colorado School of Mines (CSM)
11
+//                  United States Geological Survey (USGS)
12
+//
13
+//          Email:  tirons@mines.edu, tirons@usgs.gov
14
+//
15
+//  This program is free software: you can redistribute it and/or modify
16
+//  it under the terms of the GNU General Public License as published by
17
+//  the Free Software Foundation, either version 3 of the License, or
18
+//  (at your option) any later version.
19
+//
20
+//  This program is distributed in the hope that it will be useful,
21
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
22
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23
+//  GNU General Public License for more details.
24
+//
25
+//  You should have received a copy of the GNU General Public License
26
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
27
+//
28
+// ===========================================================================
29
+
30
+/**
31
+  @file
32
+  @author   Trevor Irons
33
+  @date     08/16/12
34
+  @version   0.0
35
+ **/
36
+
37
+#include "FEM4EllipticPDE.h"
38
+#include "vtkSphere.h"
39
+#include "vtkRectilinearGrid.h"
40
+#include "vtkFloatArray.h"
41
+#include <vtkIdList.h>
42
+#include <vtkIdTypeArray.h>
43
+#include <vtkCell.h>
44
+#include <vtkTriangleFilter.h>
45
+#include <vtkDataSetSurfaceFilter.h>
46
+#include <vtkExtractCells.h>
47
+#include <vtkGeometryFilter.h>
48
+#include <vtkUnstructuredGrid.h>
49
+#include <vtkExtractEdges.h>
50
+#include <vtkImplicitDataSet.h>
51
+#include <vtkXMLUnstructuredGridReader.h>
52
+
53
+// debugging stuff
54
+#include <vtkDataSetMapper.h>
55
+#include <vtkSelectionNode.h>
56
+#include <vtkExtractSelection.h>
57
+#include <vtkSelection.h>
58
+#include <vtkVertexGlyphFilter.h>
59
+
60
+#define RENDERTEST
61
+//#define CONNECTTEST
62
+
63
+#ifdef RENDERTEST
64
+#include "vtkRectilinearGridGeometryFilter.h"
65
+#include "vtkRectilinearGridOutlineFilter.h"
66
+#include "vtkPolyDataMapper.h"
67
+#include "vtkActor.h"
68
+#include "vtkRenderWindowInteractor.h"
69
+#include "vtkRenderer.h"
70
+#include "vtkRenderWindow.h"
71
+#include "vtkCamera.h"
72
+#include "vtkProperty.h"
73
+#include <vtkDataSetMapper.h>
74
+#include <vtkSelectionNode.h>
75
+#include <vtkSelection.h>
76
+#include <vtkExtractSelection.h>
77
+#include <vtkExtractEdges.h>
78
+#include <vtkVertexGlyphFilter.h>
79
+#include <vtkTriangleFilter.h>
80
+#include <vtkImplicitHalo.h>
81
+#endif
82
+
83
+using namespace Lemma;
84
+
85
+Real Magnet(const Real& x, const Real& y, const Real& z);
86
+
87
+int main(int argc, char**argv) {
88
+
89
+    if (argc < 4) {
90
+        std::cout << "usage: ./utFEMEllipticPDE_bhmag  <g.vtu>  <mesh.vtu>  <results.vtu>" << std::endl;
91
+        //std::cout << "usage: ./utFEMEllipticPDE_bhmag   <mesh.vtu>  <results.vtu>" << std::endl;
92
+        exit(EXIT_SUCCESS);
93
+    }
94
+
95
+    int    nx = 80;
96
+    //double dx =  .0021875 ; // 160
97
+    double dx =  .004375 ; // 160
98
+    double ox = -.175;
99
+    vtkFloatArray *xCoords = vtkFloatArray::New();
100
+    for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
101
+
102
+    int    ny = 80;
103
+    double dy =  .004375;
104
+    double oy = -.175;
105
+    vtkFloatArray *yCoords = vtkFloatArray::New();
106
+    for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
107
+
108
+    int    nz = 343; // 685; // 160
109
+    double dz = .004375;
110
+    double oz = 9.75;
111
+    vtkFloatArray *zCoords = vtkFloatArray::New();
112
+    for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
113
+
114
+    vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
115
+        rGrid->SetDimensions(nx+1, ny+1, nz+1);
116
+        rGrid->SetExtent(0, nx, 0, ny, 0,  nz);
117
+        rGrid->SetXCoordinates(xCoords);
118
+        rGrid->SetYCoordinates(yCoords);
119
+        rGrid->SetZCoordinates(zCoords);
120
+
121
+    ////////////////////////////////////////////
122
+    // Define Sigma/mu term
123
+    // TODO
124
+
125
+    ////////////////////////////////////////////
126
+    // Define source (G) term
127
+    vtkXMLUnstructuredGridReader* GReader = vtkXMLUnstructuredGridReader::New();
128
+        std::cout << "Reading G file " << argv[1] << std::endl;
129
+        GReader->SetFileName(argv[1]);
130
+        GReader->Update();
131
+
132
+    vtkImplicitDataSet* implG = vtkImplicitDataSet::New();
133
+        implG->SetDataSet(GReader->GetOutput());
134
+        implG->SetOutValue(0.);
135
+        implG->SetOutGradient(0., 0., 0.);
136
+
137
+    ////////////////////////////////////////////
138
+    // Define solution mesh
139
+    vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
140
+        std::cout << "Reading Mesh file " << argv[2] << std::endl;
141
+        MeshReader->SetFileName(argv[2]);
142
+        MeshReader->Update();
143
+
144
+    ////////////////////////////////////////////
145
+    // Solve
146
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
147
+        Solver->SetGFunction(implG);
148
+        //Solver->SetGFunction(Magnet);
149
+        //Solver->SetSigmaFunction(implSigma);
150
+        Solver->SetBoundaryStep(.05);
151
+        Solver->SetGrid(MeshReader->GetOutput());
152
+        //Solver->SetGrid(rGrid);
153
+        Solver->Solve(argv[3]);
154
+
155
+    // Clean up
156
+    Solver->Delete();
157
+    GReader->Delete();
158
+    implG->Delete();
159
+
160
+    exit(EXIT_SUCCESS);
161
+}
162
+
163
+Real Magnet(const Real& x, const Real& y, const Real& z) {
164
+
165
+    if (z < 10 || z > 11 ) {
166
+        return 0.;
167
+    } else if ( std::sqrt(x*x + y*y) <= 0.05 ) {
168
+        if (x>0.) return  1.;
169
+        if (x<0.) return -1.;
170
+    }
171
+    return 0.;
172
+}
173
+

+ 104
- 0
examples/VTKDC.cpp Näytä tiedosto

@@ -0,0 +1,104 @@
1
+/* This file is part of Lemma, a geophysical modelling and inversion API.
2
+ * More information is available at http://lemmasoftware.org
3
+ */
4
+
5
+/* This Source Code Form is subject to the terms of the Mozilla Public
6
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
7
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
+ */
9
+
10
+/**
11
+ * @file
12
+ * @date      08/07/2014 04:16:25 PM
13
+ * @version   $Id$
14
+ * @author    Trevor Irons (ti)
15
+ * @email     Trevor.Irons@xri-geo.com
16
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
17
+ * @copyright Copyright (c) 2014, Trevor Irons
18
+ */
19
+
20
+#include <vtkUnstructuredGrid.h>
21
+
22
+#include <vtkUnstructuredGridReader.h>
23
+#include <vtkUnstructuredGridWriter.h>
24
+
25
+#include <vtkGenericDataObjectReader.h>
26
+
27
+#include <vtkXMLUnstructuredGridReader.h>
28
+#include <vtkXMLUnstructuredGridWriter.h>
29
+#include <vtkDoubleArray.h>
30
+#include <vtkPointData.h>
31
+
32
+#include <cmath>
33
+#include <Eigen/Core>
34
+
35
+#include "DCSurvey.h"
36
+#include "FEM4EllipticPDE.h"
37
+
38
+#ifdef HAVE_YAMLCPP
39
+#include "yaml-cpp/yaml.h"
40
+#endif
41
+
42
+using namespace Lemma;
43
+
44
+int main(int argc, char**argv) {
45
+
46
+    std::cout << "Mesh processing routine\n";
47
+    if (argc<4) {
48
+        std::cout << "usage:\n" << "./utVTKEdge  mesh.vtu  input.yaml   results.vtu"  << std::endl;
49
+        exit(EXIT_SUCCESS);
50
+    }
51
+
52
+    vtkUnstructuredGrid* uGrid = vtkUnstructuredGrid::New();
53
+
54
+    std::string fn = argv[1];
55
+    if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
56
+        vtkGenericDataObjectReader* Reader = vtkGenericDataObjectReader::New();
57
+            Reader->SetFileName(argv[1]);
58
+            Reader->Update();
59
+        if(Reader->IsFileUnstructuredGrid()) {
60
+            std::cout << "Found Unstructured grid legacy file" << std::endl;
61
+            uGrid = Reader->GetUnstructuredGridOutput();
62
+        } else {
63
+            std::cerr << "Unknown legacy format";
64
+            exit(EXIT_FAILURE);
65
+        }
66
+
67
+    } else {
68
+        vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
69
+            std::cout << "Reading" << argv[1] << std::endl;
70
+            Reader->SetFileName(argv[1]);
71
+            Reader->Update();
72
+            uGrid = Reader->GetOutput();
73
+    }
74
+    int nc = uGrid->GetNumberOfCells()  ;
75
+    int nn = uGrid->GetNumberOfPoints()  ;
76
+
77
+    std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
78
+
79
+    std::ifstream ifstr(argv[2]);
80
+    //std::vector<YAML::Node> nodes = YAML::LoadAll(ifstr);
81
+    YAML::Node nodes = YAML::Load(ifstr);
82
+    DCSurvey* Survey = DCSurvey::DeSerialize(nodes);
83
+    //std::cout << *Survey << std::endl;
84
+
85
+
86
+    ////////////////////////////////////////////
87
+    // Solve
88
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
89
+        //Solver->SetGFunction(implG);
90
+        //Solver->SetGFunction(Magnet);
91
+        //Solver->SetSigmaFunction(implSigma);
92
+        //Solver->SetBoundaryStep(.05);
93
+
94
+        Solver->SetGrid(uGrid);
95
+        Solver->SetupDC(Survey, 0);
96
+        Solver->Solve( argv[3] );
97
+        //Solver->SetGrid(rGrid);
98
+        //Solver->Solve(argv[3]);
99
+
100
+    Survey->Delete();
101
+    Solver->Delete();
102
+
103
+    exit(EXIT_SUCCESS);
104
+}

+ 95
- 0
examples/VTKEdgeG.cpp Näytä tiedosto

@@ -0,0 +1,95 @@
1
+/* This file is part of Lemma, a geophysical modelling and inversion API.
2
+ * More information is available at http://lemmasoftware.org
3
+ */
4
+
5
+/* This Source Code Form is subject to the terms of the Mozilla Public
6
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
7
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
+ */
9
+
10
+/**
11
+ * @file
12
+ * @date      08/07/2014 04:16:25 PM
13
+ * @version   $Id$
14
+ * @author    Trevor Irons (ti)
15
+ * @email     Trevor.Irons@xri-geo.com
16
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
17
+ * @copyright Copyright (c) 2014, Trevor Irons
18
+ */
19
+
20
+#include <vtkUnstructuredGrid.h>
21
+
22
+#include <vtkUnstructuredGridReader.h>
23
+#include <vtkUnstructuredGridWriter.h>
24
+
25
+#include <vtkXMLUnstructuredGridReader.h>
26
+#include <vtkXMLUnstructuredGridWriter.h>
27
+#include <vtkDoubleArray.h>
28
+#include <vtkPointData.h>
29
+
30
+#include <cmath>
31
+
32
+
33
+int main(int argc, char**argv) {
34
+
35
+    std::cout << "Mesh processing routine\n";
36
+    if (argc<3) {
37
+        std::cout << "usage:\n" << "./utVTKEdge  inMesh.vtu outG.vtu" << std::endl;
38
+    }
39
+
40
+    vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
41
+        std::cout << "Reading" << argv[1] << std::endl;
42
+        Reader->SetFileName(argv[1]);
43
+        Reader->Update();
44
+
45
+        int nc = Reader->GetOutput()->GetNumberOfCells()  ;
46
+        int nn = Reader->GetOutput()->GetNumberOfPoints()  ;
47
+
48
+    std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
49
+
50
+    // Input magnet size TODO get from file or .geo file even?
51
+    double R = .25;          // radius of magnet
52
+    double D0 = 10;          // Top of magnet
53
+    double D1 = 11;          // Bottom of magnet
54
+    double eps = 1e-2;       // epsilon for on the point
55
+    // Pol = double[3];
56
+
57
+    // Declare new array for the data
58
+    vtkDoubleArray* G = vtkDoubleArray::New();
59
+        G->SetNumberOfComponents(1);
60
+        G->SetNumberOfTuples(nn);
61
+        G->SetName("G");
62
+
63
+    // Loop over nodes, look at position, set with G if on boundary
64
+    double point[3];
65
+    for (int in=0; in<nn; ++in) {
66
+        Reader->GetOutput()->GetPoint(in, &point[0]);
67
+        //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
68
+        double raddist = sqrt( point[0]*point[0] + point[1]*point[1] );
69
+        if ( raddist > R - eps && raddist < R + eps && point[2] < D1 && point[2] > D0 ) {
70
+            if ( point[1] < 0 ) {
71
+                G->InsertTuple1( in,   1 );
72
+            } else if (point[1] > 0) {
73
+                G->InsertTuple1( in,  -1 );
74
+            } else {
75
+                G->InsertTuple1( in, 0 );
76
+            }
77
+        } else {
78
+            G->InsertTuple1( in, 0 );
79
+        }
80
+
81
+    }
82
+
83
+    // Add new data
84
+    Reader->GetOutput()->GetPointData()->SetScalars( G );
85
+
86
+    // Write out new file
87
+    vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
88
+        Writer->SetInputConnection(Reader->GetOutputPort());
89
+        Writer->SetFileName( argv[2] );
90
+        Writer->Write();
91
+
92
+    Reader->Delete();
93
+
94
+
95
+}

+ 172
- 0
examples/VTKEdgeGsphere.cpp Näytä tiedosto

@@ -0,0 +1,172 @@
1
+/* This file is part of Lemma, a geophysical modelling and inversion API.
2
+ * More information is available at http://lemmasoftware.org
3
+ */
4
+
5
+/* This Source Code Form is subject to the terms of the Mozilla Public
6
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
7
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
+ */
9
+
10
+/**
11
+ * @file
12
+ * @date      08/07/2014 04:16:25 PM
13
+ * @version   $Id$
14
+ * @author    Trevor Irons (ti)
15
+ * @email     Trevor.Irons@xri-geo.com
16
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
17
+ * @copyright Copyright (c) 2014, Trevor Irons
18
+ */
19
+
20
+#include <vtkUnstructuredGrid.h>
21
+
22
+#include <vtkUnstructuredGridReader.h>
23
+#include <vtkUnstructuredGridWriter.h>
24
+
25
+#include <vtkGenericDataObjectReader.h>
26
+
27
+#include <vtkXMLUnstructuredGridReader.h>
28
+#include <vtkXMLUnstructuredGridWriter.h>
29
+#include <vtkDoubleArray.h>
30
+#include <vtkPointData.h>
31
+
32
+#include <cmath>
33
+#include <Eigen/Core>
34
+
35
+const double PI = 4.0*atan(1.0);
36
+
37
+int main(int argc, char**argv) {
38
+
39
+    std::cout << "Mesh processing routine\n";
40
+    if (argc<4) {
41
+        std::cout << "usage:\n" << "./utVTKEdge  inMesh.vtu  <radius> <problemType>  outG.vtu" << std::endl;
42
+        exit(EXIT_SUCCESS);
43
+    }
44
+
45
+    vtkUnstructuredGrid* uGrid = vtkUnstructuredGrid::New();
46
+
47
+    std::string fn = argv[1];
48
+    if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
49
+        vtkGenericDataObjectReader* Reader = vtkGenericDataObjectReader::New();
50
+            Reader->SetFileName(argv[1]);
51
+            Reader->Update();
52
+        if(Reader->IsFileUnstructuredGrid()) {
53
+            std::cout << "Found Unstructured grid legacy file" << std::endl;
54
+            uGrid = Reader->GetUnstructuredGridOutput();
55
+        } else {
56
+            std::cerr << "Unknown legacy format";
57
+            exit(EXIT_FAILURE);
58
+        }
59
+
60
+    } else {
61
+
62
+        vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
63
+            std::cout << "Reading" << argv[1] << std::endl;
64
+            Reader->SetFileName(argv[1]);
65
+            Reader->Update();
66
+            uGrid = Reader->GetOutput();
67
+    }
68
+
69
+
70
+        int nc = uGrid->GetNumberOfCells()  ;
71
+        int nn = uGrid->GetNumberOfPoints()  ;
72
+
73
+    std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
74
+
75
+    // Input magnet size TODO get from file or .geo file even?
76
+    double R =  atof(argv[2]); //.25;          // radius of magnet
77
+    double eps = 1e-4;       // epsilon for on the point
78
+    // Pol = double[3];
79
+
80
+    // Declare new array for the data
81
+    vtkDoubleArray* G = vtkDoubleArray::New();
82
+        G->SetNumberOfComponents(1);
83
+        G->SetNumberOfTuples(nn);
84
+        G->SetName("G");
85
+
86
+    vtkDoubleArray* phi = vtkDoubleArray::New();
87
+        phi->SetNumberOfComponents(1);
88
+        phi->SetNumberOfTuples(nn);
89
+        phi->SetName("analytic_phi");
90
+
91
+    // Loop over nodes, look at position, set with G if on boundary
92
+    double point[3];
93
+    Eigen::Vector3d M; M << 0,0,1;
94
+    for (int in=0; in<nn; ++in) {
95
+        uGrid->GetPoint(in, &point[0]);
96
+        //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
97
+        double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
98
+
99
+        //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
100
+        // Calculate RHS
101
+        if ( std::string(argv[3]) == "gravity") {
102
+            if ( raddist < R + eps) {
103
+                G->InsertTuple1( in, 1 );
104
+            } else {
105
+                G->InsertTuple1( in, 0 );
106
+            }
107
+        }
108
+
109
+        else if ( std::string(argv[3]) == "magnet") {
110
+            if ( raddist > R - eps && raddist < R + eps ) {
111
+                // \rho = \nabla \cdot \mathbf{M}
112
+                //
113
+                if ( point[2] < -eps ) {
114
+                    G->InsertTuple1( in,   1 );
115
+                } else if (point[2] > eps) {
116
+                    G->InsertTuple1( in,  -1 );
117
+                } else {
118
+                    G->InsertTuple1( in, 0 );
119
+                }
120
+                //
121
+                // Above is inaccurate and unstable
122
+
123
+                // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
124
+                Eigen::Vector3d n;
125
+                n << point[0],point[1],point[2];
126
+                n.array() /= raddist;
127
+                //double sigma = n.dot(M);
128
+                G->InsertTuple1(in, n.dot(M) );
129
+            } else {
130
+                G->InsertTuple1( in, 0 );
131
+            }
132
+        }
133
+
134
+        // Insert analytic phi part
135
+        /* magnetics problem p. 198 Jackson */
136
+        if (std::string(argv[3]) == "magnet") {
137
+            if (raddist < R) {
138
+                phi->InsertTuple1( in, (1./3.)*point[2] );
139
+            } else {
140
+                phi->InsertTuple1( in, 0);
141
+                double costheta = point[2]/raddist ;
142
+                //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
143
+                phi->InsertTuple1( in, (1./3.) * (R*R*R) * (costheta / (raddist*raddist))  );
144
+            }
145
+        } else if (std::string(argv[3]) == "gravity") {
146
+            double mass = 4./3. * PI * (R*R*R);
147
+            if (raddist < R) {
148
+                phi->InsertTuple1( in, mass * ( 3*(R*R) - raddist*raddist )/(2*(R*R*R)) ); // (1./3.)*point[2] );
149
+            } else {
150
+                //phi->InsertTuple1( in, 0);
151
+                //double costheta = point[2]/raddist ;
152
+                //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
153
+                phi->InsertTuple1( in, mass/raddist );
154
+            }
155
+        }
156
+    }
157
+
158
+    // Add new data
159
+    uGrid->GetPointData()->AddArray( phi );
160
+    uGrid->GetPointData()->SetScalars( G );
161
+
162
+    // Write out new file
163
+    vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
164
+        //Writer->SetInputConnection(Reader->GetOutputPort());
165
+        Writer->SetInputData( uGrid );
166
+        Writer->SetFileName( argv[4] );
167
+        Writer->Write();
168
+
169
+    //Reader->Delete();
170
+
171
+
172
+}

+ 27
- 0
examples/borehole/Makefile Näytä tiedosto

@@ -0,0 +1,27 @@
1
+
2
+sphereG.vtu : sphere.geo
3
+	#gmsh  -3 -format vtk -o sphere.vtk  sphere.geo
4
+	../utVTKEdgeGsphere sphere.vtk .25 gravity  sphereG.vtu
5
+	#../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
6
+
7
+sphereOut: 
8
+	../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 
9
+	#paraview --state=sliceSphere.pvsm
10
+
11
+all : sphereG.vtu sphereOut
12
+
13
+# --data=sphereOut.vtu
14
+
15
+# even finer meshing
16
+# radius = 4, u = \pm .9              analytic = \pm 1.33
17
+
18
+# Finer meshing
19
+# radius = 5, u = \pm 1.65            analytic = \pm 1.67
20
+# radius = 4, u = \pm 1.09            analytic = \pm 1.33
21
+
22
+# radius = 5, u = \pm 2.5             analytic = \pm 1.67
23
+# radius = 4, u = \pm 1.66            analytic = \pm 1.33
24
+# radius = 3, u= \pm .926 -.985       analytic = \pm 1
25
+# radius = 2, u = \pm.42              analytic = \pm .667
26
+# for radius=1,  u = \pm .111         analytic = \pm .333 
27
+# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167

+ 31
- 0
examples/borehole/Notes.txt Näytä tiedosto

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

+ 100
- 0
examples/borehole/boreholeMagnet.geo Näytä tiedosto

@@ -0,0 +1,100 @@
1
+/**
2
+ *  Trevor P. Irons, XRI Geophysics, LLC
3
+ *  Test of Coloumbic Magnetic Potential
4
+ */
5
+
6
+lc =   .01;   // Target element size
7
+R =    .25;   // Magnet Radius
8
+D0 = 10;      // Top of magnet
9
+D1 = 11;      // Bottom of magnet
10
+
11
+// Total Solution Space
12
+X0 = -1.;
13
+X1 =  1.;
14
+Y0 = -1.;
15
+Y1 =  1.;
16
+Z0 =  9.5;
17
+Z1 = 11.5;
18
+
19
+////////////////////////////////////
20
+// North Pole
21
+p = newp;
22
+dd = 5;
23
+d = lc/(1*dd);
24
+Point(p) =   { 0,  0, D0, lc/dd};
25
+Point(p+1) = { R,  d, D0, lc/dd};
26
+Point(p+2) = { R, -d, D0, lc/dd};
27
+Point(p+3) = {-R,  d, D0, lc/dd};
28
+Point(p+4) = {-R, -d, D0, lc/dd};
29
+
30
+// Connect up the points
31
+c = newc;
32
+Circle(c)  = {p+3, p, p+1};
33
+
34
+c2 = newc;
35
+Circle(c2)  = {p+4, p, p+2};
36
+
37
+ss = news;
38
+ss = Extrude {0, 0, D1-D0} { Line{c}; };
39
+
40
+ss2 = news;
41
+ss2 = Extrude {0, 0, D1-D0} { Line{c2}; };
42
+
43
+
44
+//Symmetry { 1.0, 1.0, 1.0, 0.0 }{Duplicata{Surface{ss[1]};}}
45
+//lcs = newl;
46
+//lcs = Symmetry { 0.0, 0.0,.0001, 0.0 }{Duplicata{Line{c};}};
47
+//lcs =
48
+//Symmetry { 0.0, 1.0, 0.0, 0.0 }{Duplicata{Line{c};}}
49
+
50
+//Symmetry {0, 0, 120, 0} {
51
+// Duplicata{Surface{ss[1]};}
52
+//}
53
+
54
+
55
+
56
+
57
+//ss2 = news;
58
+//ss = Extrude {0, 0, D1-D0} { Line{c, lcs}; };
59
+//ss2 = Extrude {0, 0, D1-D0} { Line{l, lcs}; };
60
+
61
+//Symmetry {1,0,0,0}{ Duplicata{Surface{ss[1]};}};
62
+//Symmetry{ expression-list }{ transform-list }:
63
+
64
+//ss2 = news;
65
+//ss2 = Extrude {0, 0, D1-D0} { Line{c2, l2}; };
66
+
67
+/////////////////////////////////////
68
+// Large Bounding box
69
+pp = newp;
70
+Point(pp)    = {X0, Y0, Z0, lc};
71
+Point(pp+1)  = {X1, Y0, Z0, lc};
72
+Point(pp+2)  = {X1, Y1, Z0, lc};
73
+Point(pp+3)  = {X0, Y1, Z0, lc};
74
+
75
+
76
+lv = newl;
77
+Line(lv) = {pp,pp+1};
78
+Line(lv+1) = {pp+1,pp+2};
79
+Line(lv+2) = {pp+2,pp+3};
80
+Line(lv+3) = {pp+3,pp};
81
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
82
+
83
+// Hard coded doom
84
+Plane Surface(125) = {lv+4};
85
+
86
+//v = newv;
87
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
88
+
89
+Surface{ss[1]} In Volume{v[1]};
90
+Surface{ss2[1]} In Volume{v[1]};
91
+
92
+Physical Volume(1) = {v[1]};
93
+
94
+//Physical Volume("minus") = {1};
95
+//Physical Volume("background") = {3};
96
+//Coherence;
97
+
98
+//Coherence;
99
+
100
+//Mesh.RecombineAll = 1;

+ 120
- 0
examples/borehole/boreholeMagnetGrid.geo Näytä tiedosto

@@ -0,0 +1,120 @@
1
+/**
2
+ *  Trevor P. Irons, XRI Geophysics, LLC
3
+ *  Test of Coloumbic Magnetic Potential
4
+ */
5
+
6
+lc = 1e-2;   // Target element size
7
+R =  .05;    // Magnet Radius
8
+D0 = 9.5;    // Top of mesh
9
+D1 = 11.5;   // Bottom of mesh
10
+
11
+// Total Solution Space
12
+X0 = -2.;
13
+X1 =  2.;
14
+Y0 = -2.;
15
+Y1 =  2.;
16
+//Z0 =  9.5;
17
+//Z1 = 11.5;
18
+
19
+////////////////////////////////////
20
+// North Pole
21
+
22
+Function Ring
23
+    X *= 1.2;
24
+    // North Pole
25
+    p0 = newp; Point(p0) = {   0,   0, D0, lc};
26
+    p1 = newp; Point(p1) = { X*R,   0, D0, lc};
27
+    p2 = newp; Point(p2) = {   0, X*R, D0, lc};
28
+    p3 = newp; Point(p3) = {-X*R,   0, D0, lc};
29
+
30
+    c1 = newc; Circle(c1) = {p1, p0, p3};
31
+    l1 = newl; Line(l1) = {p3,p1};
32
+    l2 = newl; Line Loop(l2) = {c1, l1};
33
+    s1 = news; Plane Surface(s1) = {l2};
34
+    Extrude {0, 0, D1-D0} { Surface{s1};  }
35
+
36
+    // South Pole
37
+    p4 = newp; Point(p4)  = {   0,  -0.0001, D0, lc}; //  8
38
+    p5 = newp; Point(p5)  = { X*R,  -0.0001, D0, lc}; //  9
39
+    p6 = newp; Point(p6)  = {   0,  -X*R,    D0, lc}; // 10
40
+    p7 = newp; Point(p7)  = {-X*R,  -0.0001, D0, lc}; // 11
41
+
42
+    c2 = newc; Circle(c2) = {p7, p4, p5};
43
+    l3 = newl; Line(l3) = {p5,p7};
44
+    l4 = newl; Line Loop(l4) = {c2, l3};
45
+    s2 = news; Plane Surface(s2) = {l4};
46
+    Extrude {0, 0, D1-D0} { Surface{s2};  }
47
+
48
+    //Compound Volume(4) = {vol2[1], vol3[1]};
49
+
50
+//Circle(12) = {11, 8, 9};
51
+//Line(13) = {9, 11};
52
+//Line Loop(14) = {12, 13};
53
+//Plane Surface(15) = {14};
54
+
55
+Return
56
+
57
+X = 1.;
58
+For r In {1:10}
59
+  Call Ring;
60
+EndFor
61
+
62
+/*
63
+
64
+Point(0) = { 0, 0, D0, lc};
65
+Point(1) = { R, 0, D0, lc};
66
+Point(2) = { 0, R, D0, lc};
67
+Point(3) = {-R, 0, D0, lc};
68
+
69
+// Connect up the points
70
+Circle(4) = {1, 0, 3};
71
+Line(5) = {3,1};
72
+Line Loop(6) = {4,5};
73
+Plane Surface(7) = {6};
74
+
75
+////////////////////////////////////
76
+// South Pole
77
+
78
+Point(8)  = { 0,  -0.0001, D0, lc};
79
+Point(9)  = { R,  -0.0001, D0, lc};
80
+Point(10) = { 0, -R, D0, lc};
81
+Point(11) = {-R,  -0.0001, D0, lc};
82
+
83
+// Connect up the points
84
+Circle(12) = {11, 8, 9};
85
+Line(13) = {9, 11};
86
+Line Loop(14) = {12, 13};
87
+Plane Surface(15) = {14};
88
+
89
+//////////////////////////////////////
90
+// Extrude magnet
91
+Extrude {0, 0, D1-D0} { Surface{7};  }
92
+Extrude {0, 0, D1-D0} { Surface{15}; }
93
+
94
+/////////////////////////////////////
95
+// Large Bounding box
96
+Point(116)  = {X0, Y0, Z0};
97
+Point(117)  = {X1, Y0, Z0};
98
+Point(118)  = {X1, Y1, Z0};
99
+Point(119)  = {X0, Y1, Z0};
100
+Line(120) = {116,117};
101
+Line(121) = {117,118};
102
+Line(122) = {118,119};
103
+Line(123) = {119,116};
104
+Line Loop(124) = {120, 121, 122, 123};
105
+Plane Surface(125) = {124};
106
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
107
+
108
+//////////////////////////////////////
109
+// Rings of sensitivity calculation
110
+//Circle(222) = {2.*1, 0, 2.*2};
111
+
112
+
113
+//////////////////////////////////////
114
+// Volumes
115
+
116
+Physical Volume("plus") = {2};
117
+Physical Volume("minus") = {1};
118
+Physical Volume("background") = {3};
119
+Coherence;
120
+*/

+ 85
- 0
examples/borehole/boreholeMagnetGrid_merge.geo Näytä tiedosto

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

+ 75
- 0
examples/borehole/boreholeMagnetMinus.geo Näytä tiedosto

@@ -0,0 +1,75 @@
1
+/**
2
+ *  Trevor P. Irons, XRI Geophysics, LLC
3
+ *  Test of Coloumbic Magnetic Potential
4
+ */
5
+
6
+lc = 1e-2;   // Target element size
7
+R = .05;   // Magnet Radius
8
+D0 = 10;   // Top of magnet
9
+D1 = 11;   // Bottom of magnet
10
+
11
+// Total Solution Space
12
+X0 = -2.;
13
+X1 =  2.;
14
+Y0 = -2.;
15
+Y1 =  2.;
16
+Z0 =  9.;
17
+Z1 = 12.;
18
+
19
+////////////////////////////////////
20
+// North Pole
21
+/*
22
+Point(0) = { 0, 0, D0, lc};
23
+Point(1) = { R, 0, D0, lc};
24
+Point(2) = { 0, R, D0, lc};
25
+Point(3) = {-R, 0, D0, lc};
26
+
27
+// Connect up the points
28
+Circle(4) = {1, 0, 3};
29
+Line(5) = {3,1};
30
+Line Loop(6) = {4,5};
31
+Plane Surface(7) = {6};
32
+*/
33
+
34
+
35
+////////////////////////////////////
36
+// South Pole
37
+
38
+Point(8)  = { 0,  -0.0001, D0, lc};
39
+Point(9)  = { R,  -0.0001, D0, lc};
40
+Point(10) = { 0, -R, D0, lc};
41
+Point(11) = {-R,  -0.0001, D0, lc};
42
+
43
+// Connect up the points
44
+Circle(12) = {11, 8, 9};
45
+Line(13) = {9, 11};
46
+Line Loop(14) = {12, 13};
47
+Plane Surface(15) = {14};
48
+
49
+
50
+//////////////////////////////////////
51
+// Extrude magnet
52
+//Extrude {0, 0, D1-D0} { Surface{7};  }
53
+Extrude {0, 0, D1-D0} { Surface{15}; }
54
+
55
+/*
56
+// Large Bounding box
57
+Point(116)  = {X0, Y0, Z0};
58
+Point(117)  = {X1, Y0, Z0};
59
+Point(118)  = {X1, Y1, Z0};
60
+Point(119)  = {X0, Y1, Z0};
61
+Line(120) = {116,117};
62
+Line(121) = {117,118};
63
+Line(122) = {118,119};
64
+Line(123) = {119,116};
65
+Line Loop(124) = {120, 121, 122, 123}; 
66
+Plane Surface(125) = {124};
67
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
68
+*/
69
+
70
+//////////////////////////////////////
71
+// Volumes
72
+
73
+Physical Volume("plus") = {2};
74
+Physical Volume("minus") = {1};
75
+//Physical Volume("background") = {3};

+ 175
- 0
examples/borehole/boreholeMagnetPlus.geo Näytä tiedosto

@@ -0,0 +1,175 @@
1
+/**
2
+ *  Trevor P. Irons, XRI Geophysics, LLC
3
+ *  Test of Coloumbic Magnetic Potential
4
+ */
5
+
6
+
7
+Function Cyl
8
+    // centre
9
+    pp0 = newp; Point(pp0) = {   0,    0, D0, lc2};
10
+
11
+    // Points defining outer ring
12
+    pp1 = newp; Point(pp1) = { R2,   0, D0, lc2};
13
+    pp2 = newp; Point(pp2) = {  0,  R2, D0, lc2};
14
+    pp3 = newp; Point(pp3) = {-R2,   0, D0, lc2};
15
+    pp4 = newp; Point(pp4) = {  0, -R2, D0, lc2};
16
+
17
+    cc1 = newc; Circle(cc1) = {pp1, pp0, pp2};
18
+    cc2 = newc; Circle(cc2) = {pp2, pp0, pp3};
19
+    cc3 = newc; Circle(cc3) = {pp3, pp0, pp4};
20
+    cc4 = newc; Circle(cc4) = {pp4, pp0, pp1};
21
+
22
+    ll1 = newl; Line Loop(ll1) = {cc1, cc2, cc3, cc4};
23
+
24
+    ss1 = news; Plane Surface(ss1) = {ll1};
25
+    //v1 = newv; v1 =
26
+    //newv vol2;
27
+    vol2[] = Extrude {0, 0, D1-D0} { Surface{ss1}; };
28
+
29
+Return
30
+
31
+lc = 5e-2;   // Target element size
32
+R = .05;   // Magnet Radius
33
+D0 = 10;   // Top of magnet
34
+D1 = 11;   // Bottom of magnet
35
+
36
+
37
+p0 = newp; Point(p0) = { 0, 0, D0, lc};
38
+p1 = newp; Point(p1) = { R, 0, D0, lc};
39
+p2 = newp; Point(p2) = { 0, R, D0, lc};
40
+p3 = newp; Point(p3) = {-R, 0, D0, lc};
41
+p4 = newp; Point(p4) = { 0,-R, D0, lc};
42
+
43
+c1 = newc; Circle(c1) = {p1, p0, p2};
44
+c2 = newc; Circle(c2) = {p2, p0, p3};
45
+c3 = newc; Circle(c3) = {p3, p0, p4};
46
+c4 = newc; Circle(c4) = {p4, p0, p1};
47
+
48
+
49
+ll0 = newl; Line Loop(ll0) = {1,2,3,4};
50
+Plane Surface(6) = {5};
51
+vol2[] = Extrude {0,0,D1-D0} {
52
+  Surface{6};
53
+};
54
+
55
+
56
+/*
57
+vol1[] = Extrude {0, 0, D1-D0} { Line{c1:c4}; };
58
+*/
59
+////////////////////////////////////////////////////
60
+// Big Cylinder
61
+////////////////////////////////////////////////////
62
+/*
63
+lc2 = 1e-2;    // Target element size
64
+R2 =  .175;    // Minimum Radius
65
+D0 = 9.75;     // Top of mesh
66
+D1 = 11.25;    // Bottom of mesh
67
+
68
+Call Cyl;
69
+
70
+//newv vol3;
71
+//Compound Volume(vol3) = {vol1[1], vol2[1]};
72
+
73
+
74
+Field[1] = Attractor;
75
+Field[1].NodesList = {p0,p1,p2,p3,p4};
76
+
77
+Field[2] = Threshold;
78
+Field[2].IField = 1;
79
+Field[2].LcMin = lc / 5;
80
+Field[2].LcMax = lc;
81
+Field[2].DistMin = 0.05;
82
+Field[2].DistMax = 0.5;
83
+
84
+// Use minimum of all the fields as the background field
85
+Field[3] = Min;
86
+Field[3].FieldsList = {2};
87
+Background Field = 2;
88
+
89
+// Don't extend the elements sizes from the boundary inside the domain
90
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
91
+Mesh.CharacteristicLengthFromPoints = 1 ;
92
+*/
93
+Mesh 3;
94
+
95
+/*
96
+Delete { Volume{vol1[1]}; }
97
+Delete { Volume{vol2[1]}; }
98
+
99
+Delete { Point{p0}; }
100
+Delete { Point{p1}; }
101
+Delete { Point{p2}; }
102
+Delete { Point{p3}; }
103
+Delete { Point{p4}; }
104
+
105
+Delete { Line{c1}; }
106
+Delete { Line{c2}; }
107
+Delete { Line{c3}; }
108
+Delete { Line{c4}; }
109
+*/
110
+//Delete { Line{ll0}; }
111
+//Delete { Volume{vol2}; }
112
+
113
+
114
+
115
+
116
+
117
+
118
+/*
119
+////////////////////////////////////
120
+// South Pole
121
+
122
+Point(8)  = { 0,  -0.0001, D0, lc};
123
+Point(9)  = { R,  -0.0001, D0, lc};
124
+Point(10) = { 0, -R, D0, lc};
125
+Point(11) = {-R,  -0.0001, D0, lc};
126
+
127
+// Connect up the points
128
+Circle(12) = {11, 8, 9};
129
+Line(13) = {9, 11};
130
+Line Loop(14) = {12, 13};
131
+Plane Surface(15) = {14};
132
+*/
133
+
134
+//////////////////////////////////////
135
+// Extrude magnet
136
+//Extrude {0, 0, D1-D0} { Surface{7};  }
137
+//Extrude {0, 0, D1-D0} { Surface{15}; }
138
+
139
+/*
140
+// Large Bounding box
141
+Point(116)  = {X0, Y0, Z0};
142
+Point(117)  = {X1, Y0, Z0};
143
+Point(118)  = {X1, Y1, Z0};
144
+Point(119)  = {X0, Y1, Z0};
145
+Line(120) = {116,117};
146
+Line(121) = {117,118};
147
+Line(122) = {118,119};
148
+Line(123) = {119,116};
149
+Line Loop(124) = {120, 121, 122, 123};
150
+Plane Surface(125) = {124};
151
+Extrude {0, 0, Z1-Z0} { Surface{125}; }
152
+*/
153
+
154
+//////////////////////////////////////
155
+// Volumes
156
+
157
+//Physical Volume("plus") = {2};
158
+//Physical Volume("minus") = {1};
159
+//Physical Volume("background") = {3};
160
+/*
161
+Delete {
162
+  Surface{9};
163
+}
164
+Delete {
165
+  Surface{13};
166
+}
167
+Delete {
168
+  Surface{17};
169
+}
170
+Delete {
171
+  Surface{21};
172
+}
173
+*/
174
+Compound Volume(49) = {1};
175
+Compound Volume(50) = {49, 1};

+ 146
- 0
examples/borehole/sphere.geo Näytä tiedosto

@@ -0,0 +1,146 @@
1
+/* This file is part of Lemma, a geophysical modelling and inversion API.
2
+ * More information is available at http://lemmasoftware.org
3
+ */
4
+
5
+/* This Source Code Form is subject to the terms of the Mozilla Public
6
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
7
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
+ */
9
+
10
+/**
11
+ * @file
12
+ * @date      08/08/2014 12:19:20 PM
13
+ * @version   $Id$
14
+ * @author    Trevor Irons (ti)
15
+ * @email     Trevor.Irons@xri-geo.com
16
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
17
+ * @copyright Copyright (c) 2014, Trevor Irons
18
+ */
19
+
20
+
21
+D0 = 10;          // Top of magnet
22
+D1 = 11;          // Bottom of magnet
23
+radius = 2.25;     // Radius of the damn thing
24
+
25
+lc = radius*10;   //  0.25;   // Target element size
26
+
27
+// Total Solution Space
28
+Box = 10*radius; // The down side of potential
29
+X0 = -Box;
30
+X1 =  Box;
31
+Y0 = -Box;
32
+Y1 =  Box;
33
+Z0 = -Box;
34
+Z1 =  Box;
35
+
36
+cellSize=lc; //300;
37
+dd = 0 ; //  1e-5; //cellSize; // .01;
38
+pio2=Pi/2;
39
+
40
+// Calculate offset effect
41
+theta = Asin(dd/radius);
42
+rr = radius * Cos(theta);
43
+
44
+///////////////////////////////////
45
+// Positive half sphere
46
+// create inner 1/8 shell
47
+p0 = newp;
48
+Point(p0)   = {      0,   0,       0, cellSize}; // origin
49
+Point(p0+1) = {    -rr,   0,      dd, cellSize};
50
+Point(p0+2) = {      0,  rr,      dd, cellSize};
51
+Point(p0+3) = {      0,   0,  radius, cellSize};
52
+Point(p0+4) = {      0,   0,      dd, cellSize}; // origin
53
+
54
+c0 = newc;
55
+Circle(c0  ) = {p0+1, p0+4, p0+2};       // Tricky, This one needs to be offset!
56
+Circle(c0+1) = {p0+3, p0, p0+1};
57
+Circle(c0+2) = {p0+3, p0, p0+2};
58
+
59
+Line Loop(10) = {c0, -(c0+2), c0+1} ;
60
+Ruled Surface (60) = {10};
61
+
62
+////////////////////////////////////////////////////////////
63
+// Negative half sphere
64
+p = newp;
65
+Point(  p) = {      0,      0,            0, cellSize};
66
+Point(p+1) = {    -rr,      0,          -dd, cellSize};
67
+Point(p+2) = {      0,     rr,          -dd, cellSize};
68
+Point(p+3) = {      0,      0,      -radius, cellSize};
69
+Point(p+4) = {      0,      0,          -dd, cellSize};
70
+
71
+cc = newc;
72
+Circle(cc  ) = {p+1, p+4, p+2};
73
+Circle(cc+1) = {p+3, p,   p+1};
74
+Circle(cc+2) = {p+3, p,   p+2};
75
+
76
+Circle(cc+3) = {p+3, p, p+2};
77
+Circle(cc+4) = {p+3, p, p+2};
78
+Circle(cc+5) = {p+3, p, p+2};
79
+
80
+ccl = newl;
81
+Line(ccl) = { p0+3, p+3 };
82
+
83
+Line Loop(11) = {cc, -(cc+2), cc+1} ;
84
+Ruled Surface (61) = {11};
85
+
86
+// create remaining 7/8 inner shells
87
+t1[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{60};}};
88
+t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
89
+t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
90
+//
91
+t4[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{61};}};
92
+t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
93
+t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
94
+
95
+/////////////////////////////////////
96
+// Large Bounding box
97
+pp = newp;
98
+Point(pp)    = {X0, Y0, Z0, lc};
99
+Point(pp+1)  = {X1, Y0, Z0, lc};
100
+Point(pp+2)  = {X1, Y1, Z0, lc};
101
+Point(pp+3)  = {X0, Y1, Z0, lc};
102
+
103
+
104
+lv = newl;
105
+Line(lv) = {pp,pp+1};
106
+Line(lv+1) = {pp+1,pp+2};
107
+Line(lv+2) = {pp+2,pp+3};
108
+Line(lv+3) = {pp+3,pp};
109
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
110
+
111
+// Hard coded doom
112
+Plane Surface(125) = {lv+4};
113
+
114
+//v = newv;
115
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
116
+
117
+/* This is GOOD */
118
+Surface{60} In Volume{v[1]};
119
+Surface{t1[0]} In Volume{v[1]};
120
+Surface{t2[0]} In Volume{v[1]};
121
+Surface{t3[0]} In Volume{v[1]};
122
+
123
+Surface{61} In Volume{v[1]};
124
+Surface{t4[0]} In Volume{v[1]};
125
+Surface{t5[0]} In Volume{v[1]};
126
+Surface{t6[0]} In Volume{v[1]};
127
+
128
+///////////////////////////////////////////////
129
+// Attractor Field
130
+
131
+Field[1] = Attractor;
132
+Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
133
+
134
+Field[2] = MathEval;
135
+//Field[2].F = Sprintf("(.25 - F1)^2 + %g", cellSize );  // WORKS
136
+Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, cellSize/2. );
137
+Background Field = 2;
138
+
139
+// Don't extend the elements sizes from the boundary inside the domain
140
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
141
+
142
+Physical Volume(1) = {v[1]};
143
+
144
+// To create the mesh run
145
+// gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
146
+//gmsh -3 -format msh1 -o outfile.msh sphere.geo

+ 26
- 0
examples/borehole/sphere.sh Näytä tiedosto

@@ -0,0 +1,26 @@
1
+
2
+#gmsh -3 -format vtk -o sphere.vtk  sphere.geo
3
+#../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
4
+#../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 
5
+
6
+# Gravity problem
7
+../utVTKEdgeGsphere sphere.vtk .25 gravity  sphereGrav.vtu
8
+../utFEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
9
+#paraview --state=sliceSphere.pvsm
10
+
11
+
12
+# --data=sphereOut.vtu
13
+
14
+# even finer meshing
15
+# radius = 4, u = \pm .9              analytic = \pm 1.33
16
+
17
+# Finer meshing
18
+# radius = 5, u = \pm 1.65            analytic = \pm 1.67
19
+# radius = 4, u = \pm 1.09            analytic = \pm 1.33
20
+
21
+# radius = 5, u = \pm 2.5             analytic = \pm 1.67
22
+# radius = 4, u = \pm 1.66            analytic = \pm 1.33
23
+# radius = 3, u= \pm .926 -.985       analytic = \pm 1
24
+# radius = 2, u = \pm.42              analytic = \pm .667
25
+# for radius=1,  u = \pm .111         analytic = \pm .333 
26
+# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167

+ 1
- 0
examples/borehole/untitled.geo Näytä tiedosto

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

+ 126
- 0
examples/merge.cpp Näytä tiedosto

@@ -0,0 +1,126 @@
1
+/* This file is part of Lemma, a geophysical modelling and inversion API.
2
+ * More information is available at http://lemmasoftware.org
3
+ */
4
+
5
+/* This Source Code Form is subject to the terms of the Mozilla Public
6
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
7
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
+ */
9
+
10
+/**
11
+ * @file
12
+ * @date      07/29/2014 02:01:21 PM
13
+ * @version   $Id$
14
+ * @author    Trevor Irons (ti)
15
+ * @email     Trevor.Irons@xri-geo.com
16
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
17
+ * @copyright Copyright (c) 2014, Trevor Irons
18
+ */
19
+
20
+#include <vtkVersion.h>
21
+#include <vtkSmartPointer.h>
22
+#include <vtkPolyData.h>
23
+#include <vtkSphereSource.h>
24
+#include <vtkConeSource.h>
25
+#include <vtkXMLPolyDataReader.h>
26
+#include <vtkCleanPolyData.h>
27
+#include <vtkAppendPolyData.h>
28
+#include <vtkPolyDataMapper.h>
29
+#include <vtkActor.h>
30
+#include <vtkRenderWindow.h>
31
+#include <vtkRenderer.h>
32
+#include <vtkRenderWindowInteractor.h>
33
+
34
+int main(int argc, char *argv[])
35
+{
36
+  vtkSmartPointer<vtkPolyData> input1 =
37
+    vtkSmartPointer<vtkPolyData>::New();
38
+  vtkSmartPointer<vtkPolyData> input2 =
39
+    vtkSmartPointer<vtkPolyData>::New();
40
+
41
+  if(argc == 1) //command line arguments not specified
42
+    {
43
+    vtkSmartPointer<vtkSphereSource> sphereSource =
44
+      vtkSmartPointer<vtkSphereSource>::New();
45
+    sphereSource->SetCenter(5,0,0);
46
+    sphereSource->Update();
47
+
48
+    input1->ShallowCopy(sphereSource->GetOutput());
49
+
50
+    vtkSmartPointer<vtkConeSource> coneSource =
51
+      vtkSmartPointer<vtkConeSource>::New();
52
+    coneSource->Update();
53
+
54
+    input2->ShallowCopy(coneSource->GetOutput());
55
+    }
56
+  else
57
+  {
58
+    if(argc != 3)
59
+      {
60
+      std::cout << "argc = " << argc << std::endl;
61
+      std::cout << "Required arguments: File1 File2" << std::endl;
62
+      return EXIT_FAILURE;
63
+      }
64
+    std::string inputFilename1 = argv[1];
65
+    std::string inputFilename2 = argv[2];
66
+    vtkSmartPointer<vtkXMLPolyDataReader> reader1 =
67
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
68
+    reader1->SetFileName(inputFilename1.c_str());
69
+    reader1->Update();
70
+    input1->ShallowCopy(reader1->GetOutput());
71
+
72
+    vtkSmartPointer<vtkXMLPolyDataReader> reader2 =
73
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
74
+    reader2->SetFileName(inputFilename2.c_str());
75
+    reader2->Update();
76
+    input2->ShallowCopy(reader2->GetOutput());
77
+  }
78
+
79
+  //Append the two meshes
80
+  vtkSmartPointer<vtkAppendPolyData> appendFilter =
81
+    vtkSmartPointer<vtkAppendPolyData>::New();
82
+#if VTK_MAJOR_VERSION <= 5
83
+  appendFilter->AddInputConnection(input1->GetProducerPort());
84
+  appendFilter->AddInputConnection(input2->GetProducerPort());
85
+#else
86
+  appendFilter->AddInputData(input1);
87
+  appendFilter->AddInputData(input2);
88
+#endif
89
+  appendFilter->Update();
90
+
91
+  // Remove any duplicate points.
92
+  vtkSmartPointer<vtkCleanPolyData> cleanFilter =
93
+    vtkSmartPointer<vtkCleanPolyData>::New();
94
+  cleanFilter->SetInputConnection(appendFilter->GetOutputPort());
95
+  cleanFilter->Update();
96
+
97
+  //Create a mapper and actor
98
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
99
+    vtkSmartPointer<vtkPolyDataMapper>::New();
100
+  mapper->SetInputConnection(cleanFilter->GetOutputPort());
101
+
102
+  vtkSmartPointer<vtkActor> actor =
103
+    vtkSmartPointer<vtkActor>::New();
104
+  actor->SetMapper(mapper);
105
+
106
+  //Create a renderer, render window, and interactor
107
+  vtkSmartPointer<vtkRenderer> renderer =
108
+    vtkSmartPointer<vtkRenderer>::New();
109
+  vtkSmartPointer<vtkRenderWindow> renderWindow =
110
+    vtkSmartPointer<vtkRenderWindow>::New();
111
+  renderWindow->AddRenderer(renderer);
112
+  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
113
+    vtkSmartPointer<vtkRenderWindowInteractor>::New();
114
+  renderWindowInteractor->SetRenderWindow(renderWindow);
115
+
116
+  //Add the actors to the scene
117
+  renderer->AddActor(actor);
118
+  renderer->SetBackground(.3, .2, .1); // Background color dark red
119
+
120
+  //Render and interact
121
+  renderWindow->Render();
122
+  renderWindowInteractor->Start();
123
+
124
+  return EXIT_SUCCESS;
125
+}
126
+

Loading…
Peruuta
Tallenna