v0.14.0
bernstein_bezier_generate_base.cpp
Go to the documentation of this file.
1 /**
2  * \file bernstein_bezier_generate_base.cpp
3  * \example bernstein_bezier_generate_base.cpp
4  *
5  * Genarte and check Bernstein-Bezier base. Test validates three properties
6  * of BB polynomials from \cite ainsworth2011bernstein.
7  *
8  * 1) Multiplication
9  * 2) Integration
10  * 3) Derivative
11  *
12  */
13 
14 #include <MoFEM.hpp>
15 
16 #ifdef __cplusplus
17 extern "C" {
18 #endif
19 #include <cblas.h>
20 #include <quad.h>
21 #ifdef __cplusplus
22 }
23 #endif
24 
25 using namespace MoFEM;
26 
27 static char help[] = "...\n\n";
28 
29 int main(int argc, char *argv[]) {
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 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::Tools::diffShapeFunMBEDGE
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:68
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
EntityHandle
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
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
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::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::BernsteinBezier::generateIndicesTriTri
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:124
MoFEM.hpp
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::CoreTmp
Definition: Core.hpp:36
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM::BernsteinBezier::generateIndicesEdgeTri
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:99
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:128
QUAD_::npoints
int npoints
Definition: quad.h:29
main
int main(int argc, char *argv[])
Definition: bernstein_bezier_generate_base.cpp:29
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::BernsteinBezier::generateIndicesVertexEdge
static MoFEMErrorCode generateIndicesVertexEdge(const int N, int *alpha)
Definition: BernsteinBezier.cpp:77
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
Range
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
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:203
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
help
static char help[]
Definition: bernstein_bezier_generate_base.cpp:27
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
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
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
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_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
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::generateIndicesVertexTri
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:90
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::BernsteinBezier::generateIndicesEdgeEdge
static MoFEMErrorCode generateIndicesEdgeEdge(const int N, int *alpha)
Definition: BernsteinBezier.cpp:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
quad.h
Generated by Doxygen 1.8.17 and hosted at University of Glasgow