v0.13.2
Loading...
Searching...
No Matches
AnalyticalSurfaces.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#pragma once
15
16struct ArbitrarySurfaceData {
17 VectorDouble p;
18 VectorDouble f;
19 MatrixDouble jac;
20 bool useSuperQuadric;
21 double R1; // c - big radius
22 double rr; // small radius - thicknness
23 double tt; // t - roundness
24
25 SNES sNES;
26 KSP ksp;
27 PC pc;
28 Vec vecX, vecR;
29 Mat J;
30 SNESLineSearch lineSearch;
31
32 ArbitrarySurfaceData(const double radius_total, const double radius_small,
33 const double roundness)
34 : useSuperQuadric(true) {
35 R1 = radius_total - radius_small;
36 rr = radius_small;
37 tt = roundness;
38 jac.resize(4, 4, false);
39 f.resize(4, false);
40 p.resize(3, false);
41
42 ierr = createSNES();
43 CHKERRABORT(PETSC_COMM_WORLD, ierr);
44 }
45
46 ArbitrarySurfaceData(const double radius_total, const double radius_small)
47 : useSuperQuadric(false) {
48 R1 = radius_total - radius_small;
49 rr = radius_small;
50 tt = 1;
51 jac.resize(4, 4, false);
52 f.resize(4, false);
53 p.resize(3, false);
54 // ierr = createSNES();
55 // CHKERRABORT(PETSC_COMM_WORLD, ierr);
56 }
57
58 ArbitrarySurfaceData() = delete;
59
60 ~ArbitrarySurfaceData() {
61 if (useSuperQuadric)
62 ierr = destroySNES();
63 CHKERRABORT(PETSC_COMM_WORLD, ierr);
64 }
65
66 template <typename T1>
67 inline Tensor1<double, 3> getClosestPointToToroid(Tensor1<T1, 3> &my_point) {
68
69 VectorDouble output;
70 VectorDouble pp({my_point(0), my_point(1), my_point(2)});
71 // CHKERR solveTorusClosestPoint(pp, output);
72 if (this->useSuperQuadric) // FIXME: make it template
73 CHKERR solveTorusClosestPointSNES(pp, output);
74 else
75 CHKERR solveTorusClosestPoint(pp, output);
76
77 Tensor1<double, 3> closest_pt(output(0), output(1), output(2));
78
79 return closest_pt;
80 };
81
82 template <typename T1>
83 inline Tensor1<double, 3> getClosestPointToCone(Tensor1<T1, 3> &my_point) {
84
85 VectorDouble output;
86 VectorDouble pp({my_point(0), my_point(1), my_point(2)});
87 CHKERR solveConeClosestPoint(pp, output);
88
89 Tensor1<double, 3> closest_pt(output(0), output(1), output(2));
90
91 return closest_pt;
92 };
93
94 template <typename T1>
95 inline Tensor1<double, 3>
96 getClosestPointToSuperquadric(Tensor1<T1, 3> &my_point) {
97
98 VectorDouble output;
99 VectorDouble pp({my_point(0), my_point(1), my_point(2)});
100 CHKERR solveTorusClosestPointSNES(pp, output);
101 Tensor1<double, 3> closest_pt(output(0), output(1), output(2));
102
103 return closest_pt;
104 };
105
106 template <typename T> inline bool isPointInsideTorus(const T &var) {
107 const double &x = var(0);
108 const double &y = var(1);
109 const double &z = var(2);
110 return -pow(rr, 2) + pow(-R1 + sqrt(pow(x, 2) + pow(y, 2)), 2. / tt) +
111 pow(z, 2. / tt) <
112 0;
113 };
114
115 template <typename T> inline bool isPointInsideCone(const T &var) {
116 const double &x = var(0);
117 const double &y = var(1);
118 const double &z = var(2);
119 const double prod1 = sqrt(pow(x, 2) + pow(y, 2));
120 return -prod1 + rr * z > 0;
121 };
122
123 MoFEMErrorCode computeTorusF(const VectorDouble &p, const VectorDouble &var) {
124 MoFEMFunctionBegin;
125 const double &x = var(0);
126 const double &y = var(1);
127 const double &z = var(2);
128 const double &lambda = var(3);
129
130 const double prod1 = sqrt(pow(x, 2) + pow(y, 2));
131 f(0) = p(0) - x + (2 * lambda * x * (-R1 + prod1)) / prod1;
132 f(1) = p(1) - y + (2 * lambda * y * (-R1 + prod1)) / prod1;
133 f(2) = p(2) - z + 2 * lambda * z;
134 f(3) = -pow(rr, 2) + pow(-R1 + prod1, 2) + pow(z, 2);
135
136 MoFEMFunctionReturn(0);
137 }
138
139 MoFEMErrorCode computeTorusJacobian(const VectorDouble &var) {
140 MoFEMFunctionBegin;
141 const double &x = var(0);
142 const double &y = var(1);
143 const double &z = var(2);
144 const double &lambda = var(3);
145 const double prod0 = pow(x, 2) + pow(y, 2);
146 const double prod1 = sqrt(prod0);
147
148 jac(0, 0) = -1 + (2 * lambda * pow(x, 2)) / (prod0) -
149 (2 * lambda * pow(x, 2) * (-R1 + prod1)) / pow(prod0, 1.5) +
150 (2 * lambda * (-R1 + prod1)) / prod1;
151 jac(0, 1) = jac(1, 0) =
152 (2 * lambda * x * y) / (prod0) -
153 (2 * lambda * x * y * (-R1 + prod1)) / pow(prod0, 1.5);
154 jac(0, 2) = jac(2, 0) = 0;
155 jac(0, 3) = jac(3, 0) = (2 * x * (-R1 + prod1)) / prod1;
156 jac(1, 1) = -1 + (2 * lambda * pow(y, 2)) / (prod0) -
157 (2 * lambda * pow(y, 2) * (-R1 + prod1)) / pow(prod0, 1.5) +
158 (2 * lambda * (-R1 + prod1)) / prod1;
159 jac(1, 2) = jac(2, 1) = 0;
160 jac(1, 3) = jac(3, 1) = (2 * y * (-R1 + prod1)) / prod1;
161 jac(2, 2) = -1 + 2 * lambda;
162 jac(2, 3) = jac(3, 2) = 2 * z;
163 jac(3, 3) = 0;
164
165 MoFEMFunctionReturn(0);
166 }
167
168 MoFEMErrorCode computeConeF(const VectorDouble &p, const VectorDouble &var) {
169 MoFEMFunctionBegin;
170 const double &x = var(0);
171 const double &y = var(1);
172 const double &z = var(2);
173 const double &lambda = var(3);
174
175 const double prod1 = sqrt(pow(x, 2) + pow(y, 2));
176 f(0) = p(0) - x + (lambda * x) / prod1;
177 f(1) = p(1) - y + (lambda * y) / prod1;
178 f(2) = p(2) - lambda * rr - z;
179 f(3) = prod1 - rr * z;
180
181 MoFEMFunctionReturn(0);
182 }
183
184 MoFEMErrorCode computeConeJacobian(const VectorDouble &var) {
185 MoFEMFunctionBegin;
186 const double &x = var(0);
187 const double &y = var(1);
188 const double &z = var(2);
189 const double &lambda = var(3);
190 const double prod0 = pow(x, 2) + pow(y, 2);
191
192 jac(0, 0) =
193 -1 - (lambda * pow(x, 2)) / pow(prod0, 1.5) + lambda / sqrt(prod0);
194 jac(0, 1) = jac(1, 0) = -((lambda * x * y) / pow(prod0, 1.5));
195 jac(0, 2) = jac(2, 0) = 0;
196 jac(0, 3) = jac(3, 0) = x / sqrt(prod0);
197 jac(1, 1) =
198 -1 - (lambda * pow(y, 2)) / pow(prod0, 1.5) + lambda / sqrt(prod0);
199 jac(1, 2) = jac(2, 1) = 0;
200 jac(1, 3) = jac(3, 1) = y / sqrt(prod0);
201 jac(2, 2) = -1;
202 jac(2, 3) = jac(3, 2) = -rr;
203 jac(3, 3) = 0;
204
205 MoFEMFunctionReturn(0);
206 }
207
208 MoFEMErrorCode computeSuperquadricF(const VectorDouble &p,
209 const VectorDouble &var) {
210 MoFEMFunctionBegin;
211 const double &x = var(0);
212 const double &y = var(1);
213 const double &z = var(2);
214 const double &lambda = var(3);
215 const double prod1 = sqrt(pow(x, 2) + pow(y, 2));
216 const double tt2 = 2. / tt;
217
218 f(0) = p(0) - x +
219 (2 * lambda * x * (-R1 + prod1) * pow(abs(-R1 + prod1), -2 + tt2)) /
220 (tt * prod1);
221 f(1) = p(1) - y +
222 (2 * lambda * y * (-R1 + prod1) * pow(abs(-R1 + prod1), -2 + tt2)) /
223 (tt * prod1);
224 f(2) = p(2) - z + (2 * lambda * z * pow(abs(z), -2 + tt2)) / tt;
225 f(3) = -pow(rr, tt2) + pow(abs(-R1 + prod1), tt2) + pow(abs(z), tt2);
226
227 MoFEMFunctionReturn(0);
228 }
229
230 MoFEMErrorCode computeSuperquadricJacobian(const VectorDouble &var) {
231 MoFEMFunctionBegin;
232 const double &x = var(0);
233 const double &y = var(1);
234 double z = var(2);
235 const double &lambda = var(3);
236 const double tt2 = 2. / tt;
237 const double prod0 = pow(x, 2) + pow(y, 2);
238 const double prod1 = sqrt(prod0);
239 const double prod2 = abs(-R1 + prod1);
240 const double prod3 = pow(-R1 + prod1, 2);
241 const double prod4 = pow(prod2, -2 + tt2);
242
243
244 jac(0, 0) =
245 -1 +
246 (2 * lambda * (-2 + tt2) * pow(x, 2) * prod3 * pow(prod2, -4 + tt2)) /
247 (tt * (prod0)) +
248 (2 * lambda * pow(x, 2) * prod4) / (tt * (prod0)) -
249 (2 * lambda * pow(x, 2) * (-R1 + prod1) * prod4) /
250 (tt * pow(prod0, 1.5)) +
251 (2 * lambda * (-R1 + prod1) * prod4) / (tt * prod1);
252 jac(0, 1) = jac(1, 0) =
253 (2 * lambda * (-2 + tt2) * x * y * prod3 * pow(prod2, -4 + tt2)) /
254 (tt * (prod0)) +
255 (2 * lambda * x * y * prod4) / (tt * (prod0)) -
256 (2 * lambda * x * y * (-R1 + prod1) * prod4) / (tt * pow(prod0, 1.5));
257 jac(0, 2) = jac(2, 0) = 0;
258 jac(0, 3) = jac(3, 0) = (2 * x * (-R1 + prod1) * prod4) / (tt * prod1);
259 jac(1, 1) =
260 -1 +
261 (2 * lambda * (-2 + tt2) * pow(y, 2) * prod3 * pow(prod2, -4 + tt2)) /
262 (tt * (prod0)) +
263 (2 * lambda * pow(y, 2) * prod4) / (tt * (prod0)) -
264 (2 * lambda * pow(y, 2) * (-R1 + prod1) * prod4) /
265 (tt * pow(prod0, 1.5)) +
266 (2 * lambda * (-R1 + prod1) * prod4) / (tt * prod1);
267 jac(1, 2) = jac(2, 1) = 0;
268 jac(1, 3) = jac(3, 1) = (2 * y * (-R1 + prod1) * prod4) / (tt * prod1);
269 jac(2, 2) =
270 -1 +
271 (2 * lambda * (-2 + tt2) * pow(z, 2) * pow(abs(z), -4 + tt2)) / tt +
272 (2 * lambda * pow(abs(z), -2 + tt2)) / tt;
273 jac(2, 3) = jac(3, 2) = (2 * z * pow(abs(z), -2 + tt2)) / tt;
274 jac(3, 3) = 0;
275
276 MoFEMFunctionReturn(0);
277 }
278
279 MoFEMErrorCode PrintTorusDebug() {
280 MoFEMFunctionBegin;
281
282 // std::random_device rd; // obtain a random number from hardware
283 // std::mt19937 gen(rd()); // seed the generator
284 // std::uniform_real_distribution<> distr(-R1 - rr * 2,
285 // R1 + rr * 2); // define the range
286
287 auto my_rand_mn = [](double M, double N) {
288 return M + (rand() / (RAND_MAX / (N - M)));
289 };
290
291 std::ofstream afile("torus_pts.csv", std::ios::ate);
292 if (afile.is_open()) {
293 afile << "X,Y,Z\n";
294 VectorDouble p({0, 0, 0});
295 // tt = 2; // FIXME: TODO: testing
296 for (int n = 0; n != 1e6; n++) {
297 double x = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
298 double y = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
299 double z = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
300 VectorDouble coords({x, y, z});
301 if (useSuperQuadric)
302 CHKERR computeSuperquadricF(p, coords);
303 else
304 CHKERR computeTorusF(p, coords);
305
306 if (f(3) <= 0)
307 afile << x << "," << y << "," << z << "\n";
308 }
309 afile.close();
310 }
311 MoFEMFunctionReturn(0);
312 }
313
314 MoFEMErrorCode PrintTorusAnimDebug() {
315 MoFEMFunctionBegin;
316
317 auto my_rand_mn = [](double M, double N) {
318 return M + (rand() / (RAND_MAX / (N - M)));
319 };
320
321 std::ofstream afile0("torus_pts0.csv", std::ios::ate);
322 std::ofstream afile1("torus_pts1.csv", std::ios::ate);
323 if (afile0.is_open() && afile1.is_open()) {
324 afile0 << "X,Y,Z\n";
325 afile1 << "X,Y,Z\n";
326
327 for (int n = 0; n != 1e4; n++) {
328 double x = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
329 double y = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
330 double z = my_rand_mn(-R1 - rr * 2, R1 + rr * 2);
331 VectorDouble coords({x, y, z});
332
333 VectorDouble output;
334 if (useSuperQuadric)
335 CHKERR solveTorusClosestPointSNES(coords, output);
336 else
337 CHKERR solveTorusClosestPoint(coords, output);
338
339 afile0 << x << "," << y << "," << z << "\n";
340 afile1 << output(0) << "," << output(1) << "," << output(2) << "\n";
341 }
342 afile0.close();
343 afile1.close();
344 }
345 MoFEMFunctionReturn(0);
346 }
347
348 MoFEMErrorCode solveLinearSystem(MatrixDouble &jac, VectorDouble &f) {
349 MoFEMFunctionBeginHot;
350
351 int info = 0;
352 const int size = 4;
353
354 // array<int, 4> ipiv;
355 // FIXME: why is this not
356 // working????? TODO: investigate
357
358 // info = lapack_dgesv(size, size, jac.data().data(), size, ipiv.data(),
359 // f.data().data(), size);
360
361 auto inverse = [](double *A, int N) {
362 int *ipv = new int[N];
363 int lwork = N * N;
364 double *work = new double[lwork];
365 int info;
366 info = lapack_dgetrf(N, N, A, N, ipv);
367 info = lapack_dgetri(N, A, N, ipv, work, lwork);
368
369 delete[] ipv;
370 delete[] work;
371 };
372
373 inverse(jac.data().data(), 4);
374 VectorDouble f_cp(4);
375 f_cp = f;
376 noalias(f) = prod(jac, f_cp);
377
378 if (info != 0) {
379 SETERRQ1(PETSC_COMM_SELF, 1, "error solve dgesv info = %d", info);
380 }
381
382 MoFEMFunctionReturnHot(0);
383 }
384
385 MoFEMErrorCode solveTorusClosestPoint(const VectorDouble &p,
386 VectorDouble &output) {
387 MoFEMFunctionBegin;
388 VectorDouble var({p(0), p(1), p(2)});
389 var.resize(4);
390 var(3) = 1;
391 auto my_rand_mn = [](double M, double N) {
392 return M + (rand() / (RAND_MAX / (N - M)));
393 };
394 output.resize(4);
395 output.clear();
396 // var(0) = rand();
397 // var(1) = rand();
398 // var(2) = rand();
399 // var(3) = rand();
400 // var /= norm_2(var);
401 auto &x_k = output;
402 const double eps = 1e-8;
403 const int max_it = 35;
404
405 int i = 0;
406 double res_norm = 0;
407 for (; i != max_it; ++i) {
408 CHKERR computeTorusF(p, var);
409 CHKERR computeTorusJacobian(var);
410 CHKERR solveLinearSystem(jac, f);
411 noalias(x_k) = var - f;
412
413 res_norm = norm_2(x_k - var);
414 if (res_norm < eps) {
415 VectorDouble out({x_k(0), x_k(1), x_k(2)});
416 const double dist = norm_2(out - p);
417
418 break;
419 }
420 noalias(var) = x_k;
421 }
422 if (i > 30 && res_norm > eps)
423 cout << "WE HAVE A PROBLEM " << res_norm << endl;
424 MoFEMFunctionReturn(0);
425 }
426
427 MoFEMErrorCode solveConeClosestPoint(const VectorDouble &p,
428 VectorDouble &output) {
429 MoFEMFunctionBegin;
430 VectorDouble var({p(0), p(1), p(2)});
431 var.resize(4);
432 var(3) = 1;
433
434 output.resize(4);
435 output.clear();
436 auto &x_k = output;
437 const double eps = 1e-8;
438 const int max_it = 35;
439 int i = 0;
440 double res_norm = 0;
441 for (; i != max_it; ++i) {
442 CHKERR computeConeF(p, var);
443 CHKERR computeConeJacobian(var);
444 CHKERR solveLinearSystem(jac, f);
445 noalias(x_k) = var - f;
446
447 res_norm = norm_2(x_k - var);
448 if (res_norm < eps) {
449 VectorDouble out({x_k(0), x_k(1), x_k(2)});
450 const double dist = norm_2(out - p);
451
452 break;
453 }
454 noalias(var) = x_k;
455 }
456 if (i > 30 && res_norm > eps)
457 cout << "WE HAVE A PROBLEM " << res_norm << endl;
458 MoFEMFunctionReturn(0);
459 }
460
461 // PetscErrorCode (ArbitrarySurfaceData::*)(SNES, Vec , Vec , void *)
462
463 static PetscErrorCode FormFunction(SNES snes, Vec xvec, Vec fvec, void *ctx) {
464 PetscErrorCode ierr;
465 const PetscScalar *xx;
466
467 array<PetscInt, 4> idx = {0, 1, 2, 3};
468 VecGetArrayRead(xvec, &xx);
469
470 ArbitrarySurfaceData *my_ctx = (ArbitrarySurfaceData *)ctx;
471 VectorDouble &pt = my_ctx->p;
472 VectorDouble var({xx[0], xx[1], xx[2], xx[3]});
473 my_ctx->computeSuperquadricF(pt, var);
474 VecRestoreArrayRead(xvec, &xx);
475
476 VecSetValues(fvec, 4, idx.data(), &*my_ctx->f.begin(), INSERT_VALUES);
477 return 0;
478 }
479
480 static PetscErrorCode FormJacobian(SNES snes, Vec x, Mat M, Mat B,
481 void *ctx) {
482 const PetscScalar *xx;
483 PetscErrorCode ierr;
484 array<PetscInt, 4> idx = {0, 1, 2, 3};
485
486 VecGetArrayRead(x, &xx);
487
488 ArbitrarySurfaceData *my_ctx = (ArbitrarySurfaceData *)ctx;
489
490 VectorDouble var({xx[0], xx[1], xx[2], xx[3]});
491 my_ctx->computeSuperquadricJacobian(var);
492 auto &my_jac = my_ctx->jac;
493 VecRestoreArrayRead(x, &xx);
494
495 MatSetValues(B, 4, idx.data(), 4, idx.data(), &*my_jac.data().begin(),
496 INSERT_VALUES);
497 MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
498 MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
499 return 0;
500 }
501
502 MoFEMErrorCode solveTorusClosestPointSNES(const VectorDouble &pt,
503 VectorDouble &output) {
504 MoFEMFunctionBegin;
505 this->p = pt;
506
507 // PetscRandom rctx;
508 // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
509 // PetscRandomSetFromOptions(rctx);
510 // PetscRandomSetInterval(rctx, 0, 1); // FIXME:
511 // VecSetRandom(x, rctx);
512 VectorDouble init_guess;
513 CHKERR solveTorusClosestPoint(p, init_guess);
514
515 auto my_rand_mn = [](double M, double N) {
516 return M + (rand() / (RAND_MAX / (N - M)));
517 };
518 double *xa;
519 CHKERR VecGetArray(vecX, &xa);
520 xa[0] = init_guess(0);
521 xa[1] = init_guess(1);
522 xa[2] = init_guess(2);
523 xa[3] = init_guess(3);
524 CHKERR VecRestoreArray(vecX, &xa);
525
526 SNESSolve(sNES, NULL, vecX);
527
528 SNESConvergedReason reason;
529 CHKERR SNESGetConvergedReason(sNES, &reason);
530 if (reason < 0) {
531 int its;
532 CHKERR SNESGetIterationNumber(sNES, &its);
533 CHKERR PetscPrintf(PETSC_COMM_SELF,
534 "** WARNING Closest point iterations failed. Number "
535 "of Newton iterations = %D\n",
536 its);
537 }
538
539 const PetscScalar *mx;
540
541 VecGetArrayRead(vecX, &mx);
542 output.resize(3);
543 output(0) = mx[0];
544 output(1) = mx[1];
545 output(2) = mx[2];
546 VecRestoreArrayRead(vecX, &mx);
547
548 MoFEMFunctionReturn(0);
549 }
550
551 MoFEMErrorCode createSNES() {
552 MoFEMFunctionBegin;
553
554 SNESCreate(PETSC_COMM_WORLD, &sNES);
555 const PetscInt n = 4;
556 // VecCreate(PETSC_COMM_WORLD, &vecX);
557 // VecSetSizes(vecX, n, n);
558 // VecSetFromOptions(vecX);
559 MatCreateSeqDense(PETSC_COMM_SELF, n, n, &jac(0, 0), &J);
560
561 // MatCreate(PETSC_COMM_WORLD, &J);
562 // MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 4, 4);
563 // MatSetFromOptions(J);
564 // MatSetUp(J);
565
566 VecCreateSeq(PETSC_COMM_SELF, n, &vecR);
567 VecCreateSeq(PETSC_COMM_SELF, n, &vecX);
568
569 SNESSetFunction(sNES, vecR, FormFunction, (void *)this);
570 SNESSetJacobian(sNES, J, J, FormJacobian, (void *)this);
571
572 SNESGetKSP(sNES, &ksp);
573 KSPGetPC(ksp, &pc);
574 KSPSetType(ksp, KSPPREONLY);
575 PCSetType(pc, PCLU);
576 // PCSetType(pc, PCNONE);
577
578 KSPSetTolerances(ksp, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT, 20);
579 // SNESSetFromOptions(sNES);
580 // KSPSetFromOptions(ksp);
581 SNESSetTolerances(sNES, 1e-8, 1e-12, 0, 25, 1000);
582 SNESGetLineSearch(sNES, &lineSearch);
583 SNESLineSearchSetType(lineSearch, SNESLINESEARCHBT);
584 // SNESLineSearchSetType(lineSearch, SNESLINESEARCHBASIC);
585 // SNESLineSearchSetType(lineSearch, SNESLINESEARCHL2);
586
587 MoFEMFunctionReturn(0);
588 }
589
590 MoFEMErrorCode destroySNES() {
591 MoFEMFunctionBegin;
592
593 VecDestroy(&vecX);
594 VecDestroy(&vecR);
595 MatDestroy(&J);
596 SNESDestroy(&sNES);
597 KSPDestroy(&ksp);
598
599 MoFEMFunctionReturn(0);
600 }
601};
M
static Index< 'M', 3 > M
Definition: BasicFeTools.hpp:96
ierr
static PetscErrorCode ierr
Definition: base_functions.c:13
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
lambda
constexpr double lambda
Definition: free_surface.cpp:79
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
lapack_dgetrf
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:157
lapack_dgetri
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:185
T
const double T
Definition: mixed_reac_diff.cpp:158
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
N
const int N
Definition: speed_test.cpp:3
ArbitrarySurfaceData::isPointInsideTorus
bool isPointInsideTorus(const T &var)
Definition: AnalyticalSurfaces.hpp:106
ArbitrarySurfaceData::ArbitrarySurfaceData
ArbitrarySurfaceData(const double radius_total, const double radius_small, const double roundness)
Definition: AnalyticalSurfaces.hpp:32
ArbitrarySurfaceData::getClosestPointToSuperquadric
Tensor1< double, 3 > getClosestPointToSuperquadric(Tensor1< T1, 3 > &my_point)
Definition: AnalyticalSurfaces.hpp:96
ArbitrarySurfaceData::FormJacobian
static PetscErrorCode FormJacobian(SNES snes, Vec x, Mat M, Mat B, void *ctx)
Definition: AnalyticalSurfaces.hpp:480
ArbitrarySurfaceData::solveTorusClosestPointSNES
MoFEMErrorCode solveTorusClosestPointSNES(const VectorDouble &pt, VectorDouble &output)
Definition: AnalyticalSurfaces.hpp:502
ArbitrarySurfaceData::computeTorusJacobian
MoFEMErrorCode computeTorusJacobian(const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:139
ArbitrarySurfaceData::PrintTorusDebug
MoFEMErrorCode PrintTorusDebug()
Definition: AnalyticalSurfaces.hpp:279
ArbitrarySurfaceData::PrintTorusAnimDebug
MoFEMErrorCode PrintTorusAnimDebug()
Definition: AnalyticalSurfaces.hpp:314
ArbitrarySurfaceData::solveConeClosestPoint
MoFEMErrorCode solveConeClosestPoint(const VectorDouble &p, VectorDouble &output)
Definition: AnalyticalSurfaces.hpp:427
ArbitrarySurfaceData::destroySNES
MoFEMErrorCode destroySNES()
Definition: AnalyticalSurfaces.hpp:590
ArbitrarySurfaceData::getClosestPointToToroid
Tensor1< double, 3 > getClosestPointToToroid(Tensor1< T1, 3 > &my_point)
Definition: AnalyticalSurfaces.hpp:67
ArbitrarySurfaceData::solveLinearSystem
MoFEMErrorCode solveLinearSystem(MatrixDouble &jac, VectorDouble &f)
Definition: AnalyticalSurfaces.hpp:348
ArbitrarySurfaceData::computeSuperquadricJacobian
MoFEMErrorCode computeSuperquadricJacobian(const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:230
ArbitrarySurfaceData::ArbitrarySurfaceData
ArbitrarySurfaceData(const double radius_total, const double radius_small)
Definition: AnalyticalSurfaces.hpp:46
ArbitrarySurfaceData::computeConeJacobian
MoFEMErrorCode computeConeJacobian(const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:184
ArbitrarySurfaceData::isPointInsideCone
bool isPointInsideCone(const T &var)
Definition: AnalyticalSurfaces.hpp:115
ArbitrarySurfaceData::FormFunction
static PetscErrorCode FormFunction(SNES snes, Vec xvec, Vec fvec, void *ctx)
Definition: AnalyticalSurfaces.hpp:463
ArbitrarySurfaceData::computeTorusF
MoFEMErrorCode computeTorusF(const VectorDouble &p, const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:123
ArbitrarySurfaceData::solveTorusClosestPoint
MoFEMErrorCode solveTorusClosestPoint(const VectorDouble &p, VectorDouble &output)
Definition: AnalyticalSurfaces.hpp:385
ArbitrarySurfaceData::ArbitrarySurfaceData
ArbitrarySurfaceData()=delete
ArbitrarySurfaceData::computeSuperquadricF
MoFEMErrorCode computeSuperquadricF(const VectorDouble &p, const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:208
ArbitrarySurfaceData::createSNES
MoFEMErrorCode createSNES()
Definition: AnalyticalSurfaces.hpp:551
ArbitrarySurfaceData::computeConeF
MoFEMErrorCode computeConeF(const VectorDouble &p, const VectorDouble &var)
Definition: AnalyticalSurfaces.hpp:168
ArbitrarySurfaceData::getClosestPointToCone
Tensor1< double, 3 > getClosestPointToCone(Tensor1< T1, 3 > &my_point)
Definition: AnalyticalSurfaces.hpp:83
Generated by Doxygen 1.9.5 and hosted at University of Glasgow