/
main.cpp
215 lines (147 loc) · 7.92 KB
/
main.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
// This code creates a structured mesh for a 3D geometry of a mechanical part (a small and a bigger cylinder connected together).
//
// It is structured as follows:
//
// 1. The 2D projected surface of the geometry is created
// 2. The projection is extruded
// 3. A deformation function is applied to the big cylinder to get the final geometry.
// 4. A mechanical simulation is performed on the geometry
//
// In order to identify the regions needed in the simulation the following physical region numbers are used:
//
// --> 1 for the whole volume
// --> 2 for the bottom surface of the small cylinder (mechanical clamp)
// --> 3 for the surface inside the big cylinder
//
// Region number -1 is used for all construction geometries (not used in the simulation).
#include "sparselizardbase.h"
using namespace mathop;
void sparselizard(void)
{
// The physical region numbers as detailed above:
int vol = 1, clamp = 2, faceincylinder = 3;
// Define the x, y and z coordinate fields:
field x("x"), y("y"), z("z");
// Define pi:
double pi = getpi();
///// CONSTRUCTION OF THE SMALL CYLINDER with inner and outer radius 7mm and 15mm:
double smallinnerradius = 0.007, smallouterradius = 0.015;
// Define the points to construct the cylinder.
// Inputs are the x, y and z point coordinates.
shape p1("point", -1, {0,0,0});
shape p2("point", -1, {smallinnerradius*cos(45*2*pi/360), smallinnerradius*sin(45*2*pi/360), 0});
shape p3("point", -1, {smallinnerradius*cos(-45*2*pi/360),smallinnerradius*sin(-45*2*pi/360),0});
shape p4("point", -1, {smallouterradius*cos(45*2*pi/360), smallouterradius*sin(45*2*pi/360), 0});
shape p5("point", -1, {smallouterradius*cos(-45*2*pi/360),smallouterradius*sin(-45*2*pi/360),0});
// Define the arcs and the lines to construct the small cylinder.
// Inputs are the first point, last point and arc center point
// as well as the number of mesh nodes in the arc.
shape a1("arc", -1, {p2,p3,p1}, 16);
shape a2("arc", -1, {p4,p5,p1}, 16);
shape a3("arc", -1, {p3,p2,p1}, 8);
shape a4("arc", -1, {p5,p4,p1}, 8);
shape l1("line", -1, {p2,p4}, 3);
shape l2("line", -1, {p3,p5}, 3);
// Define the 2 quadrangle surfaces for the small cylinder.
// Inputs are the contour lines (following the contour clockwise or counter-clockwise).
shape q1("quadrangle", -1, {a1,l1,a2,l2});
shape q2("quadrangle", -1, {a3,l2,a4,l1});
///// CONSTRUCTION OF THE BIG CYLINDER with inner and outer radius 10mm and 30mm
// around center point on x axis at x equal 100mm.
double biginnerradius = 0.01, bigouterradius = 0.03;
// Define the points to construct the cylinder.
shape p6("point", -1, {0.1,0,0});
shape p7("point", -1, {0.1-biginnerradius*cos(45*2*pi/360), biginnerradius*sin(45*2*pi/360), 0});
shape p8("point", -1, {0.1-biginnerradius*cos(-45*2*pi/360),biginnerradius*sin(-45*2*pi/360),0});
shape p9("point", -1, {0.1-bigouterradius*cos(45*2*pi/360), bigouterradius*sin(45*2*pi/360), 0});
shape p10("point", -1, {0.1-bigouterradius*cos(-45*2*pi/360),bigouterradius*sin(-45*2*pi/360),0});
// Define the arcs and the lines to construct the big cylinder.
shape a5("arc", -1, {p7,p8,p6}, 8);
shape a6("arc", -1, {p9,p10,p6}, 8);
shape a7("arc", -1, {p8,p7,p6}, 16);
shape a8("arc", -1, {p10,p9,p6}, 16);
shape l3("line", -1, {p7,p9}, 5);
shape l4("line", -1, {p8,p10}, 5);
// Define the 2 quadrangle surfaces for the big cylinder.
shape q3("quadrangle", -1, {a5,l3,a6,l4});
shape q4("quadrangle", -1, {a7,l4,a8,l3});
///// CONSTRUCTION OF THE SURFACE BETWEEN THE TWO CYLINDERS
// Define the lines linking the 2 cylinders:
shape l5("line", -1, {p4,p9}, 12);
shape l6("line", -1, {p5,p10}, 12);
// Define the surface between the cylinders:
shape q5("quadrangle", -1, {a4,l6,a6,l5});
///// EXTRUDE THE TOP AND BOTTOM CYLINDERS AROUND THE CENTRAL PART OF THE SMALL CYLINDER:
// Extrude arguments are the physical region of the extruded volume,
// the extrusion height and the numbers of node layers in the extrusion.
// The central part in the small cylinder is 10mm thick and the two above and below 5mm.
double centralthickness = 0.01, thicknessaboveandbelow = 0.005;
// Create the bottom volume part:
shape v1 = q1.extrude(vol, -thicknessaboveandbelow, 3);
shape v2 = q2.extrude(vol, -thicknessaboveandbelow, 3);
// Create the central part:
shape v3 = q1.extrude(vol, centralthickness, 6);
shape v4 = q2.extrude(vol, centralthickness, 6);
// Create the top volume part.
// Extrude it from the top face of v3 and v4.
shape v3topface = v3.getsons()[5];
shape v4topface = v4.getsons()[5];
shape v5 = v3topface.extrude(vol, thicknessaboveandbelow, 3);
shape v6 = v4topface.extrude(vol, thicknessaboveandbelow, 3);
///// EXTRUDE THE BIG CYLINDER:
shape v7 = q3.extrude(vol, centralthickness, 6);
shape v8 = q4.extrude(vol, centralthickness, 6);
///// EXTRUDE THE SURFACE BETWEEN THE 2 CYLINDERS:
shape v9 = q5.extrude(vol, centralthickness, 6);
///// DEFORM THE BIG CYLINDER ABOVE AND BELOW:
v7.deform(0,0, 100*(bigouterradius-sqrt((x-0.1)*(x-0.1)+y*y))*(z-centralthickness/2));
v8.deform(0,0, 100*(bigouterradius-sqrt((x-0.1)*(x-0.1)+y*y))*(z-centralthickness/2));
///// GET THE BOTTOM SURFACE OF THE SMALL CYLINDER AND ASSIGN PHYSICAL REGION NUMBER 2
///// GET THE SURFACE INSIDE THE BIG CYLINDER AND ASSIGN PHYSICAL REGION NUMBER 3
// You can .write the mesh at any time during the geometry construction
// phase (see below) to easily know which son you are looking for.
shape clampface1 = v1.getsons()[5];
shape clampface2 = v2.getsons()[5];
clampface1.setphysicalregion(clamp);
clampface2.setphysicalregion(clamp);
shape faceincylinder1 = v7.getsons()[1];
shape faceincylinder2 = v8.getsons()[1];
faceincylinder1.setphysicalregion(faceincylinder);
faceincylinder2.setphysicalregion(faceincylinder);
///// LOAD THE MESH
// Add here all regions needed in the finite element simulation.
mesh mymesh({v1,v2,v3,v4,v5,v6,v7,v8,v9, clampface1,clampface2, faceincylinder1,faceincylinder2});
// You can write the mesh at any time during the geometry
// construction to easily debug and validate every line of code.
mymesh.write("mymesh.msh");
///// START THE MECHANICAL SIMULATION
// Nodal shape functions 'h1' with 3 components.
// Field u is the mechanical displacement.
field u("h1xyz");
// Use interpolation order 2 on 'vol', the whole domain:
u.setorder(vol, 2);
// Clamp on surface 'clamp' (i.e. 0 valued-Dirichlet conditions):
u.setconstraint(clamp);
// E is Young's modulus. nu is Poisson's ratio.
double E = 200e9, nu = 0.3;
formulation elasticity;
// The linear elasticity formulation is classical and thus predefined:
elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), E, nu));
// Add a surface force (to the inside of the big cylinder) to tilt the structure around the x axis:
elasticity += integral(faceincylinder, array1x3(0,10*(z-centralthickness/2),0)*tf(u));
elasticity.generate();
vec solu = solve(elasticity.A(), elasticity.b());
// Transfer the data from the solution vector to the u field:
u.setdata(vol, solu);
// Write the deflection with an order 2 interpolation. Exagerate the deflection by a large factor.
(0.3e10*u).write(vol, "u.pos", 2);
// Code validation line. Can be removed.
std::cout << (abs(compz(u)).max(vol, 2)[0] < 4.24272e-12 && abs(compz(u)).max(vol, 2)[0] > 4.2427e-12);
}
int main(void)
{
SlepcInitialize(0,{},0,0);
sparselizard();
SlepcFinalize();
return 0;
}