v0.13.2
Loading...
Searching...
No Matches
AnalyticalFun.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15namespace FractureMechanics {
16
17struct AnalyticalOptions {
18 double aLpha;
19 bool optionsInitialised;
20 PetscErrorCode getOptions() {
21
22 MoFEMFunctionBeginHot;
23 ierr =
24 PetscOptionsBegin(PETSC_COMM_WORLD, "",
25 "Get analytical boundary conditions options", "none");
26 CHKERRQ(ierr);
27 ierr = PetscOptionsScalar("-analytical_alpha", "crack angle", "", aLpha,
28 &aLpha, PETSC_NULL);
29 CHKERRQ(ierr);
30 PetscBool flg;
31 double fraction = 1;
32 ierr = PetscOptionsScalar(
33 "-analytical_alpha_fraction_pi",
34 "crack angle given as fraction of Pi aLpha = M_PI/fraction", "",
35 fraction, &fraction, &flg);
36 CHKERRQ(ierr);
37 if (flg == PETSC_TRUE) {
38 aLpha = M_PI / fraction;
39 }
40 ierr = PetscOptionsEnd();
41 CHKERRQ(ierr);
42 optionsInitialised = true;
43 MoFEMFunctionReturnHot(0);
44 }
45 AnalyticalOptions() : aLpha(0), optionsInitialised(false) {}
46};
47
48#ifndef __ANALITICAL_TRACTION__
49#define __ANALITICAL_TRACTION__
50
51struct AnalyticalForces : public NeumannForcesSurface::MethodForAnalyticalForce,
52 public AnalyticalOptions {
53
54 PetscErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords,
55 const VectorDouble3 &normal, VectorDouble3 &forces) {
56
57 MoFEMFunctionBeginHot;
58
59 // initailise options
60 if (!optionsInitialised) {
61 ierr = getOptions();
62 CHKERRQ(ierr);
63 }
64
65 const double x = coords[0];
66 const double y = coords[1];
67 const double z = coords[2];
68 double &fx = forces[0];
69 double &fy = forces[1];
70 double &fz = forces[2];
71
72 const double n0 = normal[0] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
73 pow(normal[2], 2));
74 const double n1 = normal[1] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
75 pow(normal[2], 2));
76 const double n2 = normal[2] / sqrt(pow(normal[0], 2) + pow(normal[1], 2) +
77 pow(normal[2], 2));
78
79 const double pi = M_PI;
80 const double alpha = aLpha;
81 const double theta = atan2(z, x * sin(alpha) + y * cos(alpha));
82
83 const double rxy =
84 sqrt(pow(x * sin(alpha) + y * cos(alpha), 2) + pow(z, 2));
85 const double sref1 = 1 / sqrt(2 * pi * rxy) * cos(theta / 2.) *
86 (1. - sin(theta / 2.) * sin(3. * theta / 2.));
87 const double sref2 = 1 / sqrt(2 * pi * rxy) * cos(theta / 2.) *
88 (1. + sin(theta / 2.) * sin(3. * theta / 2.));
89 const double sref3 = 1 / sqrt(2 * pi * rxy) * sin(theta / 2.) *
90 cos(theta / 2.) * cos(3. * theta / 2.);
91
92 const double scal1 = sref1 * pow(sin(alpha), 2);
93 const double scal2 = sref1 * pow(cos(alpha), 2);
94 const double scal3 = sref2;
95 const double scal4 = sref1 * sin(alpha) * cos(alpha);
96 const double scal5 = sref3 * sin(alpha);
97 const double scal6 = sref3 * cos(alpha);
98
99 // This is just example
100 fx = scal1 * n0 + scal4 * n1 + scal5 * n2;
101 fy = scal4 * n0 + scal2 * n1 + scal6 * n2;
102 fz = scal5 * n0 + scal6 * n1 + scal3 * n2;
103
104 fx *= -1;
105 fy *= -1;
106 fz *= -1;
107
108 MoFEMFunctionReturnHot(0);
109 }
110};
111
112#endif //__ANALITICAL_TRACTION__
113
114#ifndef __ANALITICAL_DISPLACEMENT__
115#define __ANALITICAL_DISPLACEMENT__
116
117struct AnalyticalDisp : public AnalyticalOptions {
118
119 vector<VectorDouble> dIsp;
120
121 virtual vector<VectorDouble> &operator()(const double X, const double Y,
122 const double Z) {
123
124 // initailise options
125 if (!optionsInitialised) {
126 getOptions();
127 }
128
129 dIsp.resize(1);
130 dIsp[0].resize(3);
131
132 double &x = dIsp[0][0];
133 double &y = dIsp[0][1];
134 double &z = dIsp[0][2];
135
136 // Current positions is given by
137
138 const double E = 1.0e5;
139 const double nu = 0;
140 const double pi = M_PI;
141 const double alpha = aLpha;
142
143 const double theta = atan2(Z, X * sin(alpha) + Y * cos(alpha));
144
145 const double ux =
146 (1 + nu) / E *
147 sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
148 (2 * pi)) *
149 cos(theta / 2.) * (3 - 4 * nu - cos(theta)) *
150 sin(alpha); // simply ux = f(Z) = Z*1e-3
151 const double uy =
152 (1 + nu) / E *
153 sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
154 (2 * pi)) *
155 cos(theta / 2.) * (3 - 4 * nu - cos(theta)) * cos(alpha);
156 const double uz =
157 (1 + nu) / E *
158 sqrt(sqrt(pow(X * sin(alpha) + Y * cos(alpha), 2) + pow(Z, 2)) /
159 (2 * pi)) *
160 sin(theta / 2.) * (3 - 4 * nu - cos(theta));
161
162 x = X + ux;
163 y = Y + uy;
164 z = Z + uz;
165
166 return dIsp;
167 }
168
169};
170
171// struct AnalyticalForces: public NeumannForcesSurface::MethodForAnalyticalForce {
172//
173// PetscErrorCode getForce(
174// const EntityHandle ent,
175// const VectorDouble3 &coords,
176// const VectorDouble3 &normal,
177// VectorDouble3 &forces
178// ) {
179// MoFEMFunctionBeginHot;
180//
181// const double x = coords[0];
182// const double y = coords[1];
183// const double z = coords[2];
184// double &fx = forces[0];
185// double &fy = forces[1];
186// double &fz = forces[2];
187//
188// // This is just example
189// fx = z*1000;
190// fy = 0;
191// fz = 0;
192//
193// // cerr << forces << endl;
194//
195// MoFEMFunctionReturnHot(0);
196// }
197//
198// };
199//
200// #endif //__ANALITICAL_TRACTION__
201//
202// #ifndef __ANALITICAL_DISPLACEMENT__
203// #define __ANALITICAL_DISPLACEMENT__
204//
205// struct AnalyticalDisp {
206//
207// vector<VectorDouble > dIsp;
208//
209// virtual vector<VectorDouble >& operator()(
210// const double X, const double Y, const double Z
211// ) {
212//
213// dIsp.resize(1);
214// dIsp[0].resize(3);
215//
216// double &x = dIsp[0][0];
217// double &y = dIsp[0][1];
218// double &z = dIsp[0][2];
219//
220// // Current positions is given by
221//
222// const double ux = Z*1e-3; // simply ux = f(Z) = Z*1e-3
223// const double uy = 0;
224// const double uz = 0;
225//
226// x = X+ux;
227// y = Y+uy;
228// z = Z+uz;
229//
230// return dIsp;
231//
232// }
233//
234// };
235
236}
237
238#endif //__ANALITICAL_DISPLACEMENT__
ierr
static PetscErrorCode ierr
Definition: base_functions.c:13
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FractureMechanics::AnalyticalDisp::operator()
virtual vector< VectorDouble > & operator()(const double X, const double Y, const double Z)
Definition: AnalyticalFun.hpp:121
FractureMechanics::AnalyticalForces::getForce
PetscErrorCode getForce(const EntityHandle ent, const VectorDouble3 &coords, const VectorDouble3 &normal, VectorDouble3 &forces)
Definition: AnalyticalFun.hpp:54
Generated by Doxygen 1.9.5 and hosted at University of Glasgow