v0.13.2
Loading...
Searching...
No Matches
Functions | Variables
bernstein_bezier_generate_base.cpp File Reference
#include <MoFEM.hpp>
#include <cblas.h>
#include <quad.h>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 29 of file bernstein_bezier_generate_base.cpp.

29 {
30
31 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
32
33 try {
34
35 moab::Core mb_instance;
36 moab::Interface &moab = mb_instance;
37
38 // Create MoFEM instance
39 MoFEM::Core core(moab);
40
41 constexpr int N = 5;
42
43 // Edge
44
45 MatrixInt edge_alpha(2 + NBEDGE_H1(N), 2);
46 CHKERR BernsteinBezier::generateIndicesVertexEdge(N, &edge_alpha(0, 0));
47 CHKERR BernsteinBezier::generateIndicesEdgeEdge(N, &edge_alpha(2, 0));
48 std::cout << "edge alpha " << edge_alpha << std::endl;
49
50 std::array<double, 6> edge_x_k = {0, 0, 0, 1, 0, 0};
51 MatrixDouble edge_x_alpha(edge_alpha.size1(), 3);
52 CHKERR BernsteinBezier::domainPoints3d(N, 2, edge_alpha.size1(),
53 &edge_alpha(0, 0), edge_x_k.data(),
54 &edge_x_alpha(0, 0));
55 std::cout << "domain points " << edge_x_alpha << endl;
56
57 const int M = 50;
58 MatrixDouble edge_base(M, edge_alpha.size1());
59 MatrixDouble edge_diff_base(M, edge_alpha.size1());
60
61 auto calc_lambda_on_edge = [](int M) {
62 MatrixDouble edge_lambda(M, 2);
63 for (size_t i = 0; i != M; ++i) {
64 double x = static_cast<double>(i) / (M - 1);
65 edge_lambda(i, 0) = N_MBEDGE0(x);
66 edge_lambda(i, 1) = N_MBEDGE1(x);
67 }
68 return edge_lambda;
69 };
70 auto edge_lambda = calc_lambda_on_edge(M);
71
72 CHKERR BernsteinBezier::baseFunctionsEdge(
73 N, M, edge_alpha.size1(), &edge_alpha(0, 0), &edge_lambda(0, 0),
74 Tools::diffShapeFunMBEDGE.data(), &edge_base(0, 0),
75 &edge_diff_base(0, 0));
76
77 auto print_edge_base = [](auto M, auto &edge_base, auto &edge_diff_base) {
78 MoFEMFunctionBegin;
79 for (size_t i = 0; i != M; ++i) {
80 double x = static_cast<double>(i) / (M - 1);
81 std::cout << "edge " << x << " ";
82 for (size_t j = 0; j != edge_base.size2(); ++j)
83 std::cout << edge_base(i, j) << " ";
84 std::cout << endl;
85 }
86
87 for (size_t i = 0; i != M; ++i) {
88 double x = static_cast<double>(i) / (M - 1);
89 std::cout << "diff_edge " << x << " ";
90 for (size_t j = 0; j != edge_diff_base.size2(); ++j)
91 std::cout << edge_diff_base(i, j) << " ";
92 std::cout << endl;
93 }
94 MoFEMFunctionReturn(0);
95 };
96 CHKERR print_edge_base(M, edge_base, edge_diff_base);
97
98 auto diag_n = [](int n) { return n * (n + 1) / 2; };
99
100 auto binomial_alpha_beta = [](int dim, const int *alpha, const int *beta) {
101 double f = boost::math::binomial_coefficient<double>(alpha[0] + beta[0],
102 alpha[0]);
103 for (int d = 1; d != dim; ++d) {
104 f *= boost::math::binomial_coefficient<double>(alpha[d] + beta[d],
105 alpha[d]);
106 }
107 return f;
108 };
109
110 auto check_property_one =
111 [N, diag_n, binomial_alpha_beta](
112 const int M, const auto &edge_alpha, const auto &edge_lambda,
113 auto &edge_base,
114 boost::function<MoFEMErrorCode(
115 const int N, const int gdim, const int n_alpha,
116 const int *alpha, const double *lambda,
117 const double *grad_lambda, double *base, double *grad_base)>
118 calc_base,
119 bool debug) {
120 MoFEMFunctionBegin;
121
122 MatrixInt edge_alpha2(diag_n(edge_alpha.size1()), edge_alpha.size2());
123 int k = 0;
124 for (int i = 0; i != edge_alpha.size1(); ++i) {
125 for (int j = i; j != edge_alpha.size1(); ++j, ++k) {
126 for (int d = 0; d != edge_alpha.size2(); ++d) {
127 edge_alpha2(k, d) = edge_alpha(i, d) + edge_alpha(j, d);
128 }
129 }
130 }
131
132 MatrixDouble edge_base2(edge_base.size1(), edge_alpha2.size1());
133 CHKERR calc_base(N + N, edge_base.size1(), edge_alpha2.size1(),
134 &edge_alpha2(0, 0), &edge_lambda(0, 0),
135 Tools::diffShapeFunMBEDGE.data(), &edge_base2(0, 0),
136 nullptr);
137
138 const double f0 = boost::math::binomial_coefficient<double>(N + N, N);
139 for (int g = 0; g != edge_base.size1(); ++g) {
140 int k = 0;
141 for (size_t i = 0; i != edge_alpha.size1(); ++i) {
142 for (size_t j = i; j != edge_alpha.size1(); ++j, ++k) {
143 const double f = binomial_alpha_beta(
144 edge_alpha.size2(), &edge_alpha(i, 0), &edge_alpha(j, 0));
145 const double b = f / f0;
146 const double B_check = b * edge_base2(g, k);
147 const double B_mult = edge_base(g, i) * edge_base(g, j);
148 const double error = std::abs(B_check - B_mult);
149 if (debug)
150 std::cout << "( " << k << " " << i << " " << j << " )"
151 << B_check << " " << B_mult << " " << error
152 << std::endl;
153 constexpr double eps = 1e-12;
154 if (error > eps)
155 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
156 "Property one is not working");
157 }
158 }
159 }
160 MoFEMFunctionReturn(0);
161 };
162
163 CHKERR check_property_one(M, edge_alpha, edge_lambda, edge_base,
164 BernsteinBezier::baseFunctionsEdge, false);
165
166 auto check_property_two_on_edge = [N](auto &edge_alpha, bool debug) {
167 MoFEMFunctionBegin;
168 const int rule = N;
169 int nb_gauss_pts = QUAD_1D_TABLE[rule]->npoints;
170 MatrixDouble gauss_pts(2, nb_gauss_pts, false);
171 cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[rule]->points[1], 2,
172 &gauss_pts(0, 0), 1);
173 cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[rule]->weights, 1,
174 &gauss_pts(1, 0), 1);
175 MatrixDouble edge_lambda(nb_gauss_pts, 2);
176 cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[rule]->points, 1,
177 &edge_lambda(0, 0), 1);
178 MatrixDouble edge_base(nb_gauss_pts, edge_alpha.size1());
179 CHKERR BernsteinBezier::baseFunctionsEdge(
180 N, nb_gauss_pts, edge_alpha.size1(), &edge_alpha(0, 0),
181 &edge_lambda(0, 0), Tools::diffShapeFunMBEDGE.data(),
182 &edge_base(0, 0), nullptr);
183
184 VectorDouble integral(edge_alpha.size1());
185 integral.clear();
186
187 for (int g = 0; g != nb_gauss_pts; ++g) {
188 for (size_t i = 0; i != edge_alpha.size1(); ++i) {
189 integral[i] += gauss_pts(1, g) * edge_base(g, i);
190 }
191 }
192
193 const double check_integral =
194 1 / boost::math::binomial_coefficient<double>(N + 1, 1);
195
196 for (auto i : integral) {
197 const double error = std::abs(i - check_integral);
198 if (debug)
199 std::cout << "edge integral " << i << " " << check_integral << " "
200 << error << endl;
201 constexpr double eps = 1e-12;
202 if (error > eps)
203 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
204 "Property two on edge is not working");
205 }
206
207 MoFEMFunctionReturn(0);
208 };
209
210 CHKERR check_property_two_on_edge(edge_alpha, false);
211
212 auto check_property_three_on_edge_derivatives =
213 [N](const int M, const auto &edge_alpha, const auto &edge_lambda,
214 const auto &edge_diff_base, bool debug) {
215 MoFEMFunctionBegin;
216
217 MatrixInt edge_alpha_diff(2 + NBEDGE_H1(N - 1), 2);
218 CHKERR BernsteinBezier::generateIndicesVertexEdge(
219 N - 1, &edge_alpha_diff(0, 0));
220 CHKERR BernsteinBezier::generateIndicesEdgeEdge(
221 N - 1, &edge_alpha_diff(2, 0));
222
223 MatrixDouble c0(2 + NBEDGE_H1(N), 2 + NBEDGE_H1(N - 1));
224 std::array<int, 2> diff0 = {1, 0};
225 CHKERR BernsteinBezier::genrateDerivativeIndicesEdges(
226 N, edge_alpha.size1(), &edge_alpha(0, 0), diff0.data(),
227 edge_alpha_diff.size1(), &edge_alpha_diff(0, 0), &c0(0, 0));
228 MatrixDouble c1(2 + NBEDGE_H1(N), 2 + NBEDGE_H1(N - 1));
229 std::array<int, 2> diff1 = {0, 1};
230 CHKERR BernsteinBezier::genrateDerivativeIndicesEdges(
231 N, edge_alpha.size1(), &edge_alpha(0, 0), diff1.data(),
232 edge_alpha_diff.size1(), &edge_alpha_diff(0, 0), &c1(0, 0));
233
234 MatrixDouble edge_base_diff(M, edge_alpha_diff.size1());
235 CHKERR BernsteinBezier::baseFunctionsEdge(
236 N - 1, M, edge_alpha_diff.size1(), &edge_alpha_diff(0, 0),
237 &edge_lambda(0, 0), Tools::diffShapeFunMBEDGE.data(),
238 &edge_base_diff(0, 0), nullptr);
239
240 for (size_t i = 0; i != M; ++i) {
241 for (size_t j = 0; j != edge_diff_base.size2(); ++j) {
242
243 double check_diff_base = 0;
244 for (int k = 0; k != edge_alpha_diff.size1(); k++) {
245 check_diff_base += c0(j, k) * edge_base_diff(i, k) *
246 Tools::diffShapeFunMBEDGE[0];
247 check_diff_base += c1(j, k) * edge_base_diff(i, k) *
248 Tools::diffShapeFunMBEDGE[1];
249 }
250 const double error =
251 std::abs(check_diff_base - edge_diff_base(i, j));
252
253 if (debug)
254 std::cout << "edge_diff_base " << check_diff_base << " "
255 << edge_diff_base(i, j) << " " << error << std::endl;
256 constexpr double eps = 1e-12;
257 if (error > eps)
258 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
259 "Property three on edge is not working");
260 }
261
262 if (debug)
263 std::cout << endl;
264 }
265
266 MoFEMFunctionReturn(0);
267 };
268
269 CHKERR check_property_three_on_edge_derivatives(M, edge_alpha, edge_lambda,
270 edge_diff_base, false);
271
272 // Triangle
273
274 MatrixInt face_alpha(3 + 3 * NBEDGE_H1(N) + NBFACETRI_H1(N), 3);
275 CHKERR BernsteinBezier::generateIndicesVertexTri(N, &face_alpha(0, 0));
276 std::array<int, 3> face_edge_n{N, N, N};
277 std::array<int *, 3> face_edge_ptr{&face_alpha(3, 0),
278 &face_alpha(3 + NBEDGE_H1(N), 0),
279 &face_alpha(3 + 2 * NBEDGE_H1(N), 0)};
280 CHKERR BernsteinBezier::generateIndicesEdgeTri(face_edge_n.data(),
281 face_edge_ptr.data());
282 CHKERR BernsteinBezier::generateIndicesTriTri(
283 N, &face_alpha(3 + 3 * NBEDGE_H1(N), 0));
284 // std::cout << "face alpha " << face_alpha << std::endl;
285
286 std::array<double, 9> face_x_k = {0, 0, 0, 1, 0, 0, 0, 1, 0};
287 MatrixDouble face_x_alpha(face_alpha.size1(), 3);
288 CHKERR BernsteinBezier::domainPoints3d(N, 3, face_alpha.size1(),
289 &face_alpha(0, 0), face_x_k.data(),
290 &face_x_alpha(0, 0));
291 // std::cout << "domain points " << face_x_alpha << std::endl << std::endl;
292
293 auto calc_lambda_on_face = [](auto &face_x_alpha) {
294 MatrixDouble face_lambda(face_x_alpha.size1(), 3);
295 for (size_t i = 0; i != face_x_alpha.size1(); ++i) {
296 face_lambda(i, 0) = 1 - face_x_alpha(i, 0) - face_x_alpha(i, 1);
297 face_lambda(i, 1) = face_x_alpha(i, 0);
298 face_lambda(i, 2) = face_x_alpha(i, 1);
299 }
300 return face_lambda;
301 };
302 auto face_lambda = calc_lambda_on_face(face_x_alpha);
303
304 MatrixDouble face_base(face_x_alpha.size1(), face_alpha.size1());
305 MatrixDouble face_diff_base(face_x_alpha.size1(), 2 * face_alpha.size1());
306 CHKERR BernsteinBezier::baseFunctionsTri(
307 N, face_x_alpha.size1(), face_alpha.size1(), &face_alpha(0, 0),
308 &face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &face_base(0, 0),
309 &face_diff_base(0, 0));
310 // std::cout << "face base " << face_base << std::endl;
311
312 CHKERR check_property_one(face_x_alpha.size1(), face_alpha, face_lambda,
313 face_base, BernsteinBezier::baseFunctionsTri,
314 false);
315
316 auto check_property_two_on_face = [N](auto &face_alpha, bool debug) {
317 MoFEMFunctionBegin;
318 const int rule = N;
319 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
320 MatrixDouble gauss_pts(3, nb_gauss_pts, false);
321 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
322 &gauss_pts(0, 0), 1);
323 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
324 &gauss_pts(1, 0), 1);
325 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
326 &gauss_pts(2, 0), 1);
327 MatrixDouble face_lambda(nb_gauss_pts, 3);
328 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1,
329 &face_lambda(0, 0), 1);
330 MatrixDouble face_base(nb_gauss_pts, face_alpha.size1());
331 CHKERR BernsteinBezier::baseFunctionsTri(
332 N, nb_gauss_pts, face_alpha.size1(), &face_alpha(0, 0),
333 &face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &face_base(0, 0),
334 nullptr);
335
336 VectorDouble integral(face_alpha.size1());
337 integral.clear();
338
339 for (int g = 0; g != nb_gauss_pts; ++g) {
340 for (size_t i = 0; i != face_alpha.size1(); ++i) {
341 integral[i] += gauss_pts(2, g) * face_base(g, i) / 2.;
342 }
343 }
344
345 const double check_integral =
346 0.5 / boost::math::binomial_coefficient<double>(N + 2, 2);
347
348 for (auto i : integral) {
349 const double error = std::abs(i - check_integral);
350 if (debug)
351 std::cout << "edge integral " << i << " " << check_integral << " "
352 << error << endl;
353 constexpr double eps = 1e-12;
354 if (error > eps)
355 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
356 "Property two on edge is not working");
357 }
358
359 MoFEMFunctionReturn(0);
360 };
361 CHKERR check_property_two_on_face(face_alpha, false);
362
363 auto check_property_three_on_face_derivatives =
364 [N](const int M, const auto &face_alpha, const auto &face_lambda,
365 const auto &face_diff_base, bool debug) {
366 MoFEMFunctionBegin;
367
368 const int nb_edge = NBEDGE_H1(N - 1);
369 const int nb_face = NBFACETRI_H1(N - 1);
370
371 MatrixInt face_alpha_diff(3 + 3 * nb_edge + nb_face, 3);
372 CHKERR BernsteinBezier::generateIndicesVertexTri(
373 N - 1, &face_alpha_diff(0, 0));
374
375 const std::array<int, 3> edge_n{N - 1, N - 1, N - 1};
376 std::array<int *, 3> edge_ptr{&face_alpha_diff(3, 0),
377 &face_alpha_diff(3 + nb_edge, 0),
378 &face_alpha_diff(3 + 2 * nb_edge, 0)};
379
380 CHKERR BernsteinBezier::generateIndicesEdgeTri(edge_n.data(),
381 edge_ptr.data());
382 CHKERR BernsteinBezier::generateIndicesTriTri(
383 N - 1, &face_alpha_diff(3 + 3 * nb_edge, 0));
384
385 MatrixDouble c0(face_alpha.size1(), face_alpha_diff.size1());
386 MatrixDouble c1(face_alpha.size1(), face_alpha_diff.size1());
387 MatrixDouble c2(face_alpha.size1(), face_alpha_diff.size1());
388
389 std::array<int, 3> diff0 = {1, 0, 0};
390 CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
391 N, face_alpha.size1(), &face_alpha(0, 0), diff0.data(),
392 face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c0(0, 0));
393 std::array<int, 3> diff1 = {0, 1, 0};
394 CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
395 N, face_alpha.size1(), &face_alpha(0, 0), diff1.data(),
396 face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c1(0, 0));
397 std::array<int, 3> diff2 = {0, 0, 1};
398 CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
399 N, face_alpha.size1(), &face_alpha(0, 0), diff2.data(),
400 face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c2(0, 0));
401
402 MatrixDouble face_base_diff(M, face_alpha_diff.size1());
403 CHKERR BernsteinBezier::baseFunctionsTri(
404 N - 1, M, face_alpha_diff.size1(), &face_alpha_diff(0, 0),
405 &face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(),
406 &face_base_diff(0, 0), nullptr);
407
408 for (size_t i = 0; i != M; ++i) {
409 for (size_t j = 0; j != face_diff_base.size2() / 2; ++j) {
410
411 VectorDouble3 check_diff_base(2);
412 check_diff_base.clear();
413 for (int k = 0; k != face_alpha_diff.size1(); ++k) {
414 for (int d = 0; d != 2; ++d) {
415 check_diff_base[d] += c0(j, k) * face_base_diff(i, k) *
416 Tools::diffShapeFunMBTRI[2 * 0 + d];
417 check_diff_base[d] += c1(j, k) * face_base_diff(i, k) *
418 Tools::diffShapeFunMBTRI[2 * 1 + d];
419 check_diff_base[d] += c2(j, k) * face_base_diff(i, k) *
420 Tools::diffShapeFunMBTRI[2 * 2 + d];
421 }
422 }
423 auto diff_base = getVectorAdaptor(&face_diff_base(i, 2 * j), 2);
424 const double error = norm_2(check_diff_base - diff_base);
425
426 if (debug)
427 std::cout << "face_diff_base " << check_diff_base << " "
428 << diff_base << " " << error << std::endl;
429 constexpr double eps = 1e-12;
430 if (error > eps)
431 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
432 "Property three on face is not working");
433 }
434
435 if (debug)
436 std::cout << endl;
437 }
438
439 MoFEMFunctionReturn(0);
440 };
441
442 CHKERR check_property_three_on_face_derivatives(
443 face_x_alpha.size1(), face_alpha, face_lambda, face_diff_base, false);
444
445 // Tetrahedron
446
447 MatrixInt tet_alpha_vec(4, 4);
448 CHKERR BernsteinBezier::generateIndicesVertexTet(N, &tet_alpha_vec(0, 0));
449
450 const int nb_dofs_on_egde = NBEDGE_H1(N);
451 std::array<MatrixInt, 6> tet_alpha_edge{
452 MatrixInt(nb_dofs_on_egde, 4), MatrixInt(nb_dofs_on_egde, 4),
453 MatrixInt(nb_dofs_on_egde, 4), MatrixInt(nb_dofs_on_egde, 4),
454 MatrixInt(nb_dofs_on_egde, 4), MatrixInt(nb_dofs_on_egde, 4)};
455 std::array<int, 6> tet_edge_n{N, N, N, N, N, N};
456 std::array<int *, 6> tet_edge_ptr{
457 &tet_alpha_edge[0](0, 0), &tet_alpha_edge[1](0, 0),
458 &tet_alpha_edge[2](0, 0), &tet_alpha_edge[3](0, 0),
459 &tet_alpha_edge[4](0, 0), &tet_alpha_edge[5](0, 0)};
460 CHKERR BernsteinBezier::generateIndicesEdgeTet(tet_edge_n.data(),
461 tet_edge_ptr.data());
462
463 const int nb_dofs_on_tri = NBFACETRI_H1(N);
464 std::array<MatrixInt, 4> tet_alpha_face{
465 MatrixInt(nb_dofs_on_tri, 4), MatrixInt(nb_dofs_on_tri, 4),
466 MatrixInt(nb_dofs_on_tri, 4), MatrixInt(nb_dofs_on_tri, 4)};
467 std::array<int, 4> tet_face_n{N, N, N, N};
468 std::array<int *, 4> tet_face_ptr{
469 &tet_alpha_face[0](0, 0), &tet_alpha_face[1](0, 0),
470 &tet_alpha_face[2](0, 0), &tet_alpha_face[3](0, 0)};
471 CHKERR BernsteinBezier::generateIndicesTriTet(tet_face_n.data(),
472 tet_face_ptr.data());
473
474 const int nb_dofs_on_tet = NBVOLUMETET_H1(N);
475 MatrixInt tet_alpha_tet(nb_dofs_on_tet, 4);
476 CHKERR BernsteinBezier::generateIndicesTetTet(N, &tet_alpha_tet(0, 0));
477
478 std::vector<MatrixInt *> alpha_tet;
479 alpha_tet.push_back(&tet_alpha_vec);
480 for (int e = 0; e != 6; ++e)
481 alpha_tet.push_back(&tet_alpha_edge[e]);
482 for (int f = 0; f != 4; ++f)
483 alpha_tet.push_back(&tet_alpha_face[f]);
484 alpha_tet.push_back(&tet_alpha_tet);
485
486 auto create_tet_mesh = [&](auto &moab_ref, Range &tets) {
487 MoFEMFunctionBegin;
488
489 std::array<double, 12> base_coords{0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
490 EntityHandle nodes[4];
491 for (int nn = 0; nn < 4; nn++)
492 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
493 EntityHandle tet;
494 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
495
496 MoFEM::CoreTmp<-1> m_core_ref(moab_ref, PETSC_COMM_SELF, -2);
497 MoFEM::Interface &m_field_ref = m_core_ref;
498 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
499 0, 3, BitRefLevel().set(0));
500 const int max_level = 3;
501 for (int ll = 0; ll != max_level; ll++) {
502 Range edges;
503 CHKERR m_field_ref.getInterface<BitRefManager>()
504 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
505 BitRefLevel().set(), MBEDGE, edges);
506 Range tets;
507 CHKERR m_field_ref.getInterface<BitRefManager>()
508 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
509 BitRefLevel(ll).set(), MBTET, tets);
510 MeshRefinement *m_ref;
511 CHKERR m_field_ref.getInterface(m_ref);
512 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
513 edges, BitRefLevel().set(ll + 1));
514 CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
515 }
516
517 CHKERR m_field_ref.getInterface<BitRefManager>()
518 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
519 BitRefLevel().set(max_level), MBTET,
520 tets);
521
522 // Use 10 node tets to print base
523 if (1) {
524 EntityHandle meshset;
525 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
526 CHKERR moab_ref.add_entities(meshset, tets);
527 CHKERR moab_ref.convert_entities(meshset, true, false, false);
528 CHKERR moab_ref.delete_entities(&meshset, 1);
529 }
530
531 MoFEMFunctionReturn(0);
532 };
533
534 Range tets;
535 CHKERR create_tet_mesh(moab, tets);
536 Range elem_nodes;
537 CHKERR moab.get_connectivity(tets, elem_nodes, false);
538 MatrixDouble coords(elem_nodes.size(), 3);
539 CHKERR moab.get_coords(elem_nodes, &coords(0, 0));
540
541 auto calc_lambda_on_tet = [](auto &coords) {
542 MatrixDouble lambda(coords.size1(), 4);
543 for (size_t i = 0; i != coords.size1(); ++i) {
544 lambda(i, 0) = 1 - coords(i, 0) - coords(i, 1) - coords(i, 2);
545 lambda(i, 1) = coords(i, 0);
546 lambda(i, 2) = coords(i, 1);
547 lambda(i, 3) = coords(i, 2);
548 }
549 return lambda;
550 };
551 auto lambda = calc_lambda_on_tet(coords);
552 // cerr << lambda << endl;
553
554 int nn = 0;
555 for (auto alpha_ptr : alpha_tet) {
556
557 MatrixDouble base(coords.size1(), alpha_ptr->size1());
558 MatrixDouble diff_base(coords.size1(), 3 * alpha_ptr->size1());
559 CHKERR BernsteinBezier::baseFunctionsTet(
560 N, coords.size1(), alpha_ptr->size1(), &(*alpha_ptr)(0, 0),
561 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &base(0, 0),
562 &diff_base(0, 0));
563
564 CHKERR check_property_one(coords.size1(), *alpha_ptr, lambda, base,
565 BernsteinBezier::baseFunctionsTet, false);
566
567 // cerr << *alpha_ptr << endl;
568 // cerr << base << endl;
569 MatrixDouble trans_base = trans(base);
570 for (int j = 0; j != base.size2(); ++j) {
571 Tag th;
572 double def_val[] = {0};
573 CHKERR moab.tag_get_handle(
574 ("base_" + boost::lexical_cast<std::string>(nn) + "_" +
575 boost::lexical_cast<std::string>(j))
576 .c_str(),
577 1, MB_TYPE_DOUBLE, th, MB_TAG_CREAT | MB_TAG_DENSE, def_val);
578 CHKERR moab.tag_set_data(th, elem_nodes, &trans_base(j, 0));
579 }
580 ++nn;
581 }
582
583 EntityHandle meshset;
584 CHKERR moab.create_meshset(MESHSET_SET, meshset);
585 CHKERR moab.add_entities(meshset, tets);
586 CHKERR moab.write_file("bb.vtk", "VTK", "", &meshset, 1);
587 }
588 CATCH_ERRORS;
589
590 CHKERR MoFEM::Core::Finalize();
591}
M
static Index< 'M', 3 > M
Definition: BasicFeTools.hpp:96
help
static char help[]
Definition: bernstein_bezier_generate_base.cpp:27
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
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
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
dim
const int dim
Definition: elec_phys_2D.cpp:12
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
lambda
constexpr double lambda
Definition: free_surface.cpp:79
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FTensor::d
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
g
constexpr double g
Definition: shallow_wave.cpp:63
N
const int N
Definition: speed_test.cpp:3
MoFEM::BernsteinBezier::baseFunctionsTri
static MoFEMErrorCode baseFunctionsTri(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:431
MoFEM::BernsteinBezier::genrateDerivativeIndicesTri
static MoFEMErrorCode genrateDerivativeIndicesTri(const int N, const int n_alpha, const int *alpha, const int *diff, const int n_alpha_diff, const int *alpha_diff, double *c)
Definition: BernsteinBezier.cpp:244
MoFEM::BernsteinBezier::generateIndicesVertexTri
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:90
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
MoFEM::BernsteinBezier::baseFunctionsTet
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:443
MoFEM::BernsteinBezier::baseFunctionsEdge
static MoFEMErrorCode baseFunctionsEdge(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:419
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:138
MoFEM::BernsteinBezier::generateIndicesEdgeEdge
static MoFEMErrorCode generateIndicesEdgeEdge(const int N, int *alpha)
Definition: BernsteinBezier.cpp:85
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:128
MoFEM::BernsteinBezier::generateIndicesTriTri
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:124
MoFEM::BernsteinBezier::domainPoints3d
static MoFEMErrorCode domainPoints3d(const int N, const int n_x, const int n_alpha, const int *alpha, const double *x_k, double *x_alpha)
Genrate BB points in 3d.
Definition: BernsteinBezier.cpp:284
MoFEM::BernsteinBezier::generateIndicesVertexEdge
static MoFEMErrorCode generateIndicesVertexEdge(const int N, int *alpha)
Definition: BernsteinBezier.cpp:77
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:203
MoFEM::BernsteinBezier::generateIndicesEdgeTri
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:99
MoFEM::BernsteinBezier::genrateDerivativeIndicesEdges
static MoFEMErrorCode genrateDerivativeIndicesEdges(const int N, const int n_alpha, const int *alpha, const int *diff, const int n_alpha_diff, const int *alpha_diff, double *c)
Definition: BernsteinBezier.cpp:237
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::Tools::diffShapeFunMBEDGE
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:68
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
QUAD_::npoints
int npoints
Definition: quad.h:29

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 27 of file bernstein_bezier_generate_base.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow