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;
50 std::array<double, 6> edge_x_k = {0, 0, 0, 1, 0, 0};
52 CHKERR BernsteinBezier::domainPoints3d(
N, 2, edge_alpha.size1(),
53 &edge_alpha(0, 0), edge_x_k.data(),
55 std::cout <<
"domain points " << edge_x_alpha << endl;
61 auto calc_lambda_on_edge = [](
int M) {
63 for (
size_t i = 0;
i !=
M; ++
i) {
64 double x =
static_cast<double>(
i) / (
M - 1);
70 auto edge_lambda = calc_lambda_on_edge(
M);
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));
77 auto print_edge_base = [](
auto M,
auto &edge_base,
auto &edge_diff_base) {
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) <<
" ";
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) <<
" ";
96 CHKERR print_edge_base(
M, edge_base, edge_diff_base);
98 auto diag_n = [](
int n) {
return n * (
n + 1) / 2; };
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],
103 for (
int d = 1;
d != dim; ++
d) {
104 f *= boost::math::binomial_coefficient<double>(alpha[
d] + beta[
d],
110 auto check_property_one =
111 [
N, diag_n, binomial_alpha_beta](
112 const int M,
const auto &edge_alpha,
const auto &edge_lambda,
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)>
122 MatrixInt edge_alpha2(diag_n(edge_alpha.size1()), edge_alpha.size2());
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);
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),
138 const double f0 = boost::math::binomial_coefficient<double>(
N +
N,
N);
139 for (
int g = 0;
g != edge_base.size1(); ++
g) {
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);
150 std::cout <<
"( " <<
k <<
" " <<
i <<
" " <<
j <<
" )"
151 << B_check <<
" " << B_mult <<
" " << error
153 constexpr
double eps = 1e-12;
156 "Property one is not working");
163 CHKERR check_property_one(
M, edge_alpha, edge_lambda, edge_base,
164 BernsteinBezier::baseFunctionsEdge,
false);
166 auto check_property_two_on_edge = [
N](
auto &edge_alpha,
bool debug) {
171 cblas_dcopy(nb_gauss_pts, &
QUAD_1D_TABLE[rule]->points[1], 2,
172 &gauss_pts(0, 0), 1);
174 &gauss_pts(1, 0), 1);
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);
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);
193 const double check_integral =
194 1 / boost::math::binomial_coefficient<double>(
N + 1, 1);
196 for (
auto i : integral) {
197 const double error = std::abs(
i - check_integral);
199 std::cout <<
"edge integral " <<
i <<
" " << check_integral <<
" "
201 constexpr
double eps = 1e-12;
204 "Property two on edge is not working");
210 CHKERR check_property_two_on_edge(edge_alpha,
false);
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) {
218 CHKERR BernsteinBezier::generateIndicesVertexEdge(
219 N - 1, &edge_alpha_diff(0, 0));
220 CHKERR BernsteinBezier::generateIndicesEdgeEdge(
221 N - 1, &edge_alpha_diff(2, 0));
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));
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));
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);
240 for (
size_t i = 0;
i !=
M; ++
i) {
241 for (
size_t j = 0;
j != edge_diff_base.size2(); ++
j) {
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];
251 std::abs(check_diff_base - edge_diff_base(
i,
j));
254 std::cout <<
"edge_diff_base " << check_diff_base <<
" "
255 << edge_diff_base(
i,
j) <<
" " << error << std::endl;
256 constexpr
double eps = 1e-12;
259 "Property three on edge is not working");
269 CHKERR check_property_three_on_edge_derivatives(
M, edge_alpha, edge_lambda,
270 edge_diff_base,
false);
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),
280 CHKERR BernsteinBezier::generateIndicesEdgeTri(face_edge_n.data(),
281 face_edge_ptr.data());
282 CHKERR BernsteinBezier::generateIndicesTriTri(
286 std::array<double, 9> face_x_k = {0, 0, 0, 1, 0, 0, 0, 1, 0};
288 CHKERR BernsteinBezier::domainPoints3d(
N, 3, face_alpha.size1(),
289 &face_alpha(0, 0), face_x_k.data(),
290 &face_x_alpha(0, 0));
293 auto calc_lambda_on_face = [](
auto &face_x_alpha) {
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);
302 auto face_lambda = calc_lambda_on_face(face_x_alpha);
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));
312 CHKERR check_property_one(face_x_alpha.size1(), face_alpha, face_lambda,
313 face_base, BernsteinBezier::baseFunctionsTri,
316 auto check_property_two_on_face = [
N](
auto &face_alpha,
bool debug) {
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);
326 &gauss_pts(2, 0), 1);
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),
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.;
345 const double check_integral =
346 0.5 / boost::math::binomial_coefficient<double>(
N + 2, 2);
348 for (
auto i : integral) {
349 const double error = std::abs(
i - check_integral);
351 std::cout <<
"edge integral " <<
i <<
" " << check_integral <<
" "
353 constexpr
double eps = 1e-12;
356 "Property two on edge is not working");
361 CHKERR check_property_two_on_face(face_alpha,
false);
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) {
371 MatrixInt face_alpha_diff(3 + 3 * nb_edge + nb_face, 3);
372 CHKERR BernsteinBezier::generateIndicesVertexTri(
373 N - 1, &face_alpha_diff(0, 0));
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)};
380 CHKERR BernsteinBezier::generateIndicesEdgeTri(edge_n.data(),
382 CHKERR BernsteinBezier::generateIndicesTriTri(
383 N - 1, &face_alpha_diff(3 + 3 * nb_edge, 0));
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());
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));
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);
408 for (
size_t i = 0;
i !=
M; ++
i) {
409 for (
size_t j = 0;
j != face_diff_base.size2() / 2; ++
j) {
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];
424 const double error = norm_2(check_diff_base - diff_base);
427 std::cout <<
"face_diff_base " << check_diff_base <<
" "
428 << diff_base <<
" " << error << std::endl;
429 constexpr
double eps = 1e-12;
432 "Property three on face is not working");
442 CHKERR check_property_three_on_face_derivatives(
443 face_x_alpha.size1(), face_alpha, face_lambda, face_diff_base,
false);
448 CHKERR BernsteinBezier::generateIndicesVertexTet(
N, &tet_alpha_vec(0, 0));
451 std::array<MatrixInt, 6> tet_alpha_edge{
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());
464 std::array<MatrixInt, 4> tet_alpha_face{
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());
475 MatrixInt tet_alpha_tet(nb_dofs_on_tet, 4);
476 CHKERR BernsteinBezier::generateIndicesTetTet(
N, &tet_alpha_tet(0, 0));
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);
486 auto create_tet_mesh = [&](
auto &moab_ref,
Range &tets) {
489 std::array<double, 12> base_coords{0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
491 for (
int nn = 0; nn < 4; nn++)
492 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
494 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
500 const int max_level = 3;
501 for (
int ll = 0; ll != max_level; ll++) {
504 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
508 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
512 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
518 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
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);
535 CHKERR create_tet_mesh(moab, tets);
537 CHKERR moab.get_connectivity(tets, elem_nodes,
false);
539 CHKERR moab.get_coords(elem_nodes, &coords(0, 0));
541 auto calc_lambda_on_tet = [](
auto &coords) {
543 for (
size_t i = 0;
i != coords.size1(); ++
i) {
544 lambda(
i, 0) = 1 - coords(
i, 0) - coords(
i, 1) - coords(
i, 2);
551 auto lambda = calc_lambda_on_tet(coords);
555 for (
auto alpha_ptr : alpha_tet) {
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),
564 CHKERR check_property_one(coords.size1(), *alpha_ptr,
lambda, base,
565 BernsteinBezier::baseFunctionsTet,
false);
570 for (
int j = 0;
j != base.size2(); ++
j) {
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))
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));
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);