|
|
@@ -0,0 +1,148 @@
|
|
|
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
|
+radius = 3.25; // Radius of the damn thing
|
|
|
11
|
+lc = radius/5; // 0.25; // Target element size
|
|
|
12
|
+
|
|
|
13
|
+// Total Solution Space
|
|
|
14
|
+Box = 3*radius; // The down side of potential
|
|
|
15
|
+X0 = -Box;
|
|
|
16
|
+X1 = Box;
|
|
|
17
|
+Y0 = -Box;
|
|
|
18
|
+Y1 = Box;
|
|
|
19
|
+Z0 = -Box;
|
|
|
20
|
+Z1 = Box;
|
|
|
21
|
+
|
|
|
22
|
+cellSize=radius/10; ///10;
|
|
|
23
|
+dd = 0 ; // 1e-5; //cellSize; // .01;
|
|
|
24
|
+pio2=Pi/2;
|
|
|
25
|
+
|
|
|
26
|
+
|
|
|
27
|
+/////////////////////////////////////
|
|
|
28
|
+// Large Bounding box
|
|
|
29
|
+pp = newp;
|
|
|
30
|
+Point(pp) = {X0, Y0, Z0, lc};
|
|
|
31
|
+Point(pp+1) = {X1, Y0, Z0, lc};
|
|
|
32
|
+Point(pp+2) = {X1, Y1, Z0, lc};
|
|
|
33
|
+Point(pp+3) = {X0, Y1, Z0, lc};
|
|
|
34
|
+
|
|
|
35
|
+
|
|
|
36
|
+lv = newl;
|
|
|
37
|
+Line(lv) = {pp,pp+1};
|
|
|
38
|
+Line(lv+1) = {pp+1,pp+2};
|
|
|
39
|
+Line(lv+2) = {pp+2,pp+3};
|
|
|
40
|
+Line(lv+3) = {pp+3,pp};
|
|
|
41
|
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
|
|
|
42
|
+
|
|
|
43
|
+// Hard coded doom
|
|
|
44
|
+Plane Surface(125) = {lv+4};
|
|
|
45
|
+
|
|
|
46
|
+//v = newv;
|
|
|
47
|
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
|
|
|
48
|
+
|
|
|
49
|
+// Calculate offset effect
|
|
|
50
|
+theta = Asin(dd/radius);
|
|
|
51
|
+rr = radius * Cos(theta);
|
|
|
52
|
+
|
|
|
53
|
+///////////////////////////////////
|
|
|
54
|
+// Positive half sphere
|
|
|
55
|
+// create inner 1/8 shell
|
|
|
56
|
+p0 = newp;
|
|
|
57
|
+Point(p0) = { 0, 0, 0, cellSize}; // origin
|
|
|
58
|
+Point(p0+1) = { -rr, 0, dd, cellSize};
|
|
|
59
|
+Point(p0+2) = { 0, rr, dd, cellSize};
|
|
|
60
|
+Point(p0+3) = { 0, 0, radius, cellSize};
|
|
|
61
|
+Point(p0+4) = { 0, 0, dd, cellSize}; // origin
|
|
|
62
|
+
|
|
|
63
|
+c0 = newc;
|
|
|
64
|
+Circle(c0 ) = {p0+1, p0+4, p0+2}; // Tricky, This one needs to be offset!
|
|
|
65
|
+Circle(c0+1) = {p0+3, p0, p0+1};
|
|
|
66
|
+Circle(c0+2) = {p0+3, p0, p0+2};
|
|
|
67
|
+
|
|
|
68
|
+Line Loop(10) = {c0, -(c0+2), c0+1} ;
|
|
|
69
|
+Ruled Surface (60) = {10};
|
|
|
70
|
+
|
|
|
71
|
+////////////////////////////////////////////////////////////
|
|
|
72
|
+// Negative half sphere
|
|
|
73
|
+p = newp;
|
|
|
74
|
+Point( p) = { 0, 0, 0, cellSize};
|
|
|
75
|
+Point(p+1) = { -rr, 0, -dd, cellSize};
|
|
|
76
|
+Point(p+2) = { 0, rr, -dd, cellSize};
|
|
|
77
|
+Point(p+3) = { 0, 0, -radius, cellSize};
|
|
|
78
|
+Point(p+4) = { 0, 0, -dd, cellSize};
|
|
|
79
|
+
|
|
|
80
|
+cc = newc;
|
|
|
81
|
+Circle(cc ) = {p+1, p+4, p+2};
|
|
|
82
|
+Circle(cc+1) = {p+3, p, p+1};
|
|
|
83
|
+Circle(cc+2) = {p+3, p, p+2};
|
|
|
84
|
+
|
|
|
85
|
+Circle(cc+3) = {p+3, p, p+2};
|
|
|
86
|
+Circle(cc+4) = {p+3, p, p+2};
|
|
|
87
|
+Circle(cc+5) = {p+3, p, p+2};
|
|
|
88
|
+
|
|
|
89
|
+ccl = newl;
|
|
|
90
|
+Line(ccl) = { p0+3, p+3 };
|
|
|
91
|
+
|
|
|
92
|
+Line Loop(11) = {cc, -(cc+2), cc+1} ;
|
|
|
93
|
+Ruled Surface (61) = {11};
|
|
|
94
|
+
|
|
|
95
|
+// create remaining 7/8 inner shells
|
|
|
96
|
+t1[] = Rotate {{0,0,1},{0,0,0},pio2 } {Duplicata{Surface{60};}};
|
|
|
97
|
+t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
|
|
|
98
|
+t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
|
|
|
99
|
+//
|
|
|
100
|
+t4[] = Rotate {{0,0,1},{0,0,0},pio2 } {Duplicata{Surface{61};}};
|
|
|
101
|
+t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
|
|
|
102
|
+t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
|
|
|
103
|
+
|
|
|
104
|
+/* This is GOOD */
|
|
|
105
|
+Surface{60} In Volume{v[1]};
|
|
|
106
|
+Surface{t1[0]} In Volume{v[1]};
|
|
|
107
|
+Surface{t2[0]} In Volume{v[1]};
|
|
|
108
|
+Surface{t3[0]} In Volume{v[1]};
|
|
|
109
|
+
|
|
|
110
|
+Surface{61} In Volume{v[1]};
|
|
|
111
|
+Surface{t4[0]} In Volume{v[1]};
|
|
|
112
|
+Surface{t5[0]} In Volume{v[1]};
|
|
|
113
|
+Surface{t6[0]} In Volume{v[1]};
|
|
|
114
|
+
|
|
|
115
|
+///////////////////////////////////////////////
|
|
|
116
|
+// Attractor Field
|
|
|
117
|
+
|
|
|
118
|
+// Field[1] = Attractor;
|
|
|
119
|
+// Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
|
|
|
120
|
+//
|
|
|
121
|
+// Field[2] = MathEval;
|
|
|
122
|
+// Field[2].F = Sprintf("(2.25 - F1)^1.01 + %g", radius/10 ); // WORKS
|
|
|
123
|
+//
|
|
|
124
|
+// Field[3] = MathEval;
|
|
|
125
|
+// Field[3].F = Sprintf("(12.25 - F1)^1.01 + %g", radius/5 ); // WORKS
|
|
|
126
|
+//
|
|
|
127
|
+// Field[4] = MathEval;
|
|
|
128
|
+// Field[4].F = Sprintf("(22.25 - F1)^1.01 + %g", radius ); // WORKS
|
|
|
129
|
+//
|
|
|
130
|
+// Field[5] = MathEval;
|
|
|
131
|
+// Field[5].F = Sprintf("(42.25 - F1)^1.01 + %g", radius*5 ); // WORKS
|
|
|
132
|
+//
|
|
|
133
|
+// //Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, 2*cellSize );
|
|
|
134
|
+// //Background Field = 2;
|
|
|
135
|
+//
|
|
|
136
|
+// // Finally, let's use the minimum of all the fields as the background mesh field
|
|
|
137
|
+// Field[7] = Min;
|
|
|
138
|
+// Field[7].FieldsList = {2, 3, 4, 5};
|
|
|
139
|
+// Background Field = 7;
|
|
|
140
|
+
|
|
|
141
|
+// Don't extend the elements sizes from the boundary inside the domain
|
|
|
142
|
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
|
|
|
143
|
+
|
|
|
144
|
+Physical Volume(1) = {v[1]};
|
|
|
145
|
+
|
|
|
146
|
+// To create the mesh run
|
|
|
147
|
+// gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
|
|
|
148
|
+//gmsh -3 -format msh1 -o outfile.msh sphere.geo
|
Trevor Irons
8 years ago