Genarte and check Bernstein-Bezier base. Test validates three properties of BB polynomials from [3].
#ifdef __cplusplus
extern "C" {
#endif
#include <cblas.h>
#ifdef __cplusplus
}
#endif
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR BernsteinBezier::generateIndicesVertexEdge(
N, &edge_alpha(0, 0));
CHKERR BernsteinBezier::generateIndicesEdgeEdge(
N, &edge_alpha(2, 0));
std::cout << "edge alpha " << edge_alpha << std::endl;
std::array<double, 6> edge_x_k = {0, 0, 0, 1, 0, 0};
CHKERR BernsteinBezier::domainPoints3d(
N, 2, edge_alpha.size1(),
&edge_alpha(0, 0), edge_x_k.data(),
&edge_x_alpha(0, 0));
std::cout << "domain points " << edge_x_alpha << endl;
auto calc_lambda_on_edge = [](
int M) {
for (
size_t i = 0;
i !=
M; ++
i) {
double x =
static_cast<double>(
i) / (
M - 1);
}
return edge_lambda;
};
auto edge_lambda = calc_lambda_on_edge(
M);
CHKERR BernsteinBezier::baseFunctionsEdge(
N,
M, edge_alpha.size1(), &edge_alpha(0, 0), &edge_lambda(0, 0),
Tools::diffShapeFunMBEDGE.data(), &edge_base(0, 0),
&edge_diff_base(0, 0));
auto print_edge_base = [](
auto M,
auto &edge_base,
auto &edge_diff_base) {
for (
size_t i = 0;
i !=
M; ++
i) {
double x =
static_cast<double>(
i) / (
M - 1);
std::cout << "edge " << x << " ";
for (
size_t j = 0;
j != edge_base.size2(); ++
j)
std::cout << edge_base(
i,
j) <<
" ";
std::cout << endl;
}
for (
size_t i = 0;
i !=
M; ++
i) {
double x =
static_cast<double>(
i) / (
M - 1);
std::cout << "diff_edge " << x << " ";
for (
size_t j = 0;
j != edge_diff_base.size2(); ++
j)
std::cout << edge_diff_base(
i,
j) <<
" ";
std::cout << endl;
}
};
CHKERR print_edge_base(
M, edge_base, edge_diff_base);
auto diag_n = [](
int n) {
return n * (
n + 1) / 2; };
auto binomial_alpha_beta = [](
int dim,
const int *
alpha,
const int *beta) {
double f = boost::math::binomial_coefficient<double>(alpha[0] + beta[0],
alpha[0]);
for (
int d = 1;
d !=
dim; ++
d) {
f *= boost::math::binomial_coefficient<double>(alpha[d] + beta[d],
alpha[d]);
}
};
auto check_property_one =
[
N, diag_n, binomial_alpha_beta](
const int M,
const auto &edge_alpha,
const auto &edge_lambda,
auto &edge_base,
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)>
calc_base,
MatrixInt edge_alpha2(diag_n(edge_alpha.size1()), edge_alpha.size2());
for (
int i = 0;
i != edge_alpha.size1(); ++
i) {
for (
int j =
i;
j != edge_alpha.size1(); ++
j, ++
k) {
for (
int d = 0;
d != edge_alpha.size2(); ++
d) {
edge_alpha2(
k, d) = edge_alpha(
i, d) + edge_alpha(
j, d);
}
}
}
MatrixDouble edge_base2(edge_base.size1(), edge_alpha2.size1());
CHKERR calc_base(
N +
N, edge_base.size1(), edge_alpha2.size1(),
&edge_alpha2(0, 0), &edge_lambda(0, 0),
Tools::diffShapeFunMBEDGE.data(), &edge_base2(0, 0),
nullptr);
const double f0 = boost::math::binomial_coefficient<double>(
N +
N,
N);
for (
int g = 0;
g != edge_base.size1(); ++
g) {
for (
size_t i = 0;
i != edge_alpha.size1(); ++
i) {
for (
size_t j =
i;
j != edge_alpha.size1(); ++
j, ++
k) {
const double f = binomial_alpha_beta(
edge_alpha.size2(), &edge_alpha(
i, 0), &edge_alpha(
j, 0));
const double B_check =
b * edge_base2(
g,
k);
const double B_mult = edge_base(
g,
i) * edge_base(
g,
j);
const double error = std::abs(B_check - B_mult);
std::cout <<
"( " <<
k <<
" " <<
i <<
" " <<
j <<
" )"
<< B_check << " " << B_mult << " " << error
<< std::endl;
constexpr double eps = 1e-12;
"Property one is not working");
}
}
}
};
CHKERR check_property_one(
M, edge_alpha, edge_lambda, edge_base,
BernsteinBezier::baseFunctionsEdge, false);
auto check_property_two_on_edge = [
N](
auto &edge_alpha,
bool debug) {
&gauss_pts(0, 0), 1);
&gauss_pts(1, 0), 1);
&edge_lambda(0, 0), 1);
CHKERR BernsteinBezier::baseFunctionsEdge(
N, nb_gauss_pts, edge_alpha.size1(), &edge_alpha(0, 0),
&edge_lambda(0, 0), Tools::diffShapeFunMBEDGE.data(),
&edge_base(0, 0), nullptr);
integral.clear();
for (
int g = 0;
g != nb_gauss_pts; ++
g) {
for (
size_t i = 0;
i != edge_alpha.size1(); ++
i) {
integral[
i] += gauss_pts(1,
g) * edge_base(
g,
i);
}
}
const double check_integral =
1 / boost::math::binomial_coefficient<double>(
N + 1, 1);
for (
auto i : integral) {
const double error = std::abs(
i - check_integral);
std::cout <<
"edge integral " <<
i <<
" " << check_integral <<
" "
<< error << endl;
constexpr double eps = 1e-12;
"Property two on edge is not working");
}
};
CHKERR check_property_two_on_edge(edge_alpha,
false);
auto check_property_three_on_edge_derivatives =
[
N](
const int M,
const auto &edge_alpha,
const auto &edge_lambda,
const auto &edge_diff_base,
bool debug) {
CHKERR BernsteinBezier::generateIndicesVertexEdge(
N - 1, &edge_alpha_diff(0, 0));
CHKERR BernsteinBezier::generateIndicesEdgeEdge(
N - 1, &edge_alpha_diff(2, 0));
std::array<int, 2> diff0 = {1, 0};
CHKERR BernsteinBezier::genrateDerivativeIndicesEdges(
N, edge_alpha.size1(), &edge_alpha(0, 0), diff0.data(),
edge_alpha_diff.size1(), &edge_alpha_diff(0, 0), &c0(0, 0));
std::array<int, 2> diff1 = {0, 1};
CHKERR BernsteinBezier::genrateDerivativeIndicesEdges(
N, edge_alpha.size1(), &edge_alpha(0, 0), diff1.data(),
edge_alpha_diff.size1(), &edge_alpha_diff(0, 0), &c1(0, 0));
CHKERR BernsteinBezier::baseFunctionsEdge(
N - 1,
M, edge_alpha_diff.size1(), &edge_alpha_diff(0, 0),
&edge_lambda(0, 0), Tools::diffShapeFunMBEDGE.data(),
&edge_base_diff(0, 0), nullptr);
for (
size_t i = 0;
i !=
M; ++
i) {
for (
size_t j = 0;
j != edge_diff_base.size2(); ++
j) {
double check_diff_base = 0;
for (
int k = 0;
k != edge_alpha_diff.size1();
k++) {
check_diff_base += c0(
j,
k) * edge_base_diff(
i,
k) *
Tools::diffShapeFunMBEDGE[0];
check_diff_base += c1(
j,
k) * edge_base_diff(
i,
k) *
Tools::diffShapeFunMBEDGE[1];
}
const double error =
std::abs(check_diff_base - edge_diff_base(
i,
j));
std::cout << "edge_diff_base " << check_diff_base << " "
<< edge_diff_base(
i,
j) <<
" " << error << std::endl;
constexpr double eps = 1e-12;
"Property three on edge is not working");
}
std::cout << endl;
}
};
CHKERR check_property_three_on_edge_derivatives(
M, edge_alpha, edge_lambda,
edge_diff_base, false);
CHKERR BernsteinBezier::generateIndicesVertexTri(
N, &face_alpha(0, 0));
std::array<int, 3> face_edge_n{
N,
N,
N};
std::array<int *, 3> face_edge_ptr{&face_alpha(3, 0),
CHKERR BernsteinBezier::generateIndicesEdgeTri(face_edge_n.data(),
face_edge_ptr.data());
CHKERR BernsteinBezier::generateIndicesTriTri(
std::array<double, 9> face_x_k = {0, 0, 0, 1, 0, 0, 0, 1, 0};
CHKERR BernsteinBezier::domainPoints3d(
N, 3, face_alpha.size1(),
&face_alpha(0, 0), face_x_k.data(),
&face_x_alpha(0, 0));
auto calc_lambda_on_face = [](auto &face_x_alpha) {
for (
size_t i = 0;
i != face_x_alpha.size1(); ++
i) {
face_lambda(
i, 0) = 1 - face_x_alpha(
i, 0) - face_x_alpha(
i, 1);
face_lambda(
i, 1) = face_x_alpha(
i, 0);
face_lambda(
i, 2) = face_x_alpha(
i, 1);
}
return face_lambda;
};
auto face_lambda = calc_lambda_on_face(face_x_alpha);
MatrixDouble face_base(face_x_alpha.size1(), face_alpha.size1());
MatrixDouble face_diff_base(face_x_alpha.size1(), 2 * face_alpha.size1());
CHKERR BernsteinBezier::baseFunctionsTri(
N, face_x_alpha.size1(), face_alpha.size1(), &face_alpha(0, 0),
&face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &face_base(0, 0),
&face_diff_base(0, 0));
CHKERR check_property_one(face_x_alpha.size1(), face_alpha, face_lambda,
face_base, BernsteinBezier::baseFunctionsTri,
false);
auto check_property_two_on_face = [
N](
auto &face_alpha,
bool debug) {
&gauss_pts(0, 0), 1);
&gauss_pts(1, 0), 1);
&gauss_pts(2, 0), 1);
&face_lambda(0, 0), 1);
CHKERR BernsteinBezier::baseFunctionsTri(
N, nb_gauss_pts, face_alpha.size1(), &face_alpha(0, 0),
&face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &face_base(0, 0),
nullptr);
integral.clear();
for (
int g = 0;
g != nb_gauss_pts; ++
g) {
for (
size_t i = 0;
i != face_alpha.size1(); ++
i) {
integral[
i] += gauss_pts(2,
g) * face_base(
g,
i) / 2.;
}
}
const double check_integral =
0.5 / boost::math::binomial_coefficient<double>(
N + 2, 2);
for (
auto i : integral) {
const double error = std::abs(
i - check_integral);
std::cout <<
"edge integral " <<
i <<
" " << check_integral <<
" "
<< error << endl;
constexpr double eps = 1e-12;
"Property two on edge is not working");
}
};
CHKERR check_property_two_on_face(face_alpha,
false);
auto check_property_three_on_face_derivatives =
[
N](
const int M,
const auto &face_alpha,
const auto &face_lambda,
const auto &face_diff_base,
bool debug) {
MatrixInt face_alpha_diff(3 + 3 * nb_edge + nb_face, 3);
CHKERR BernsteinBezier::generateIndicesVertexTri(
N - 1, &face_alpha_diff(0, 0));
const std::array<int, 3> edge_n{
N - 1,
N - 1,
N - 1};
std::array<int *, 3> edge_ptr{&face_alpha_diff(3, 0),
&face_alpha_diff(3 + nb_edge, 0),
&face_alpha_diff(3 + 2 * nb_edge, 0)};
CHKERR BernsteinBezier::generateIndicesEdgeTri(edge_n.data(),
edge_ptr.data());
CHKERR BernsteinBezier::generateIndicesTriTri(
N - 1, &face_alpha_diff(3 + 3 * nb_edge, 0));
MatrixDouble c0(face_alpha.size1(), face_alpha_diff.size1());
MatrixDouble c1(face_alpha.size1(), face_alpha_diff.size1());
MatrixDouble c2(face_alpha.size1(), face_alpha_diff.size1());
std::array<int, 3> diff0 = {1, 0, 0};
CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
N, face_alpha.size1(), &face_alpha(0, 0), diff0.data(),
face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c0(0, 0));
std::array<int, 3> diff1 = {0, 1, 0};
CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
N, face_alpha.size1(), &face_alpha(0, 0), diff1.data(),
face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c1(0, 0));
std::array<int, 3> diff2 = {0, 0, 1};
CHKERR BernsteinBezier::genrateDerivativeIndicesTri(
N, face_alpha.size1(), &face_alpha(0, 0), diff2.data(),
face_alpha_diff.size1(), &face_alpha_diff(0, 0), &c2(0, 0));
CHKERR BernsteinBezier::baseFunctionsTri(
N - 1,
M, face_alpha_diff.size1(), &face_alpha_diff(0, 0),
&face_lambda(0, 0), Tools::diffShapeFunMBTRI.data(),
&face_base_diff(0, 0), nullptr);
for (
size_t i = 0;
i !=
M; ++
i) {
for (
size_t j = 0;
j != face_diff_base.size2() / 2; ++
j) {
check_diff_base.clear();
for (
int k = 0;
k != face_alpha_diff.size1(); ++
k) {
for (
int d = 0;
d != 2; ++
d) {
check_diff_base[
d] += c0(
j,
k) * face_base_diff(
i,
k) *
Tools::diffShapeFunMBTRI[2 * 0 +
d];
check_diff_base[
d] += c1(
j,
k) * face_base_diff(
i,
k) *
Tools::diffShapeFunMBTRI[2 * 1 +
d];
check_diff_base[
d] += c2(
j,
k) * face_base_diff(
i,
k) *
Tools::diffShapeFunMBTRI[2 * 2 +
d];
}
}
const double error = norm_2(check_diff_base - diff_base);
std::cout << "face_diff_base " << check_diff_base << " "
<< diff_base << " " << error << std::endl;
constexpr double eps = 1e-12;
"Property three on face is not working");
}
std::cout << endl;
}
};
CHKERR check_property_three_on_face_derivatives(
face_x_alpha.size1(), face_alpha, face_lambda, face_diff_base, false);
CHKERR BernsteinBezier::generateIndicesVertexTet(
N, &tet_alpha_vec(0, 0));
std::array<MatrixInt, 6> tet_alpha_edge{
std::array<int, 6> tet_edge_n{
N,
N,
N,
N,
N,
N};
std::array<int *, 6> tet_edge_ptr{
&tet_alpha_edge[0](0, 0), &tet_alpha_edge[1](0, 0),
&tet_alpha_edge[2](0, 0), &tet_alpha_edge[3](0, 0),
&tet_alpha_edge[4](0, 0), &tet_alpha_edge[5](0, 0)};
CHKERR BernsteinBezier::generateIndicesEdgeTet(tet_edge_n.data(),
tet_edge_ptr.data());
std::array<MatrixInt, 4> tet_alpha_face{
std::array<int, 4> tet_face_n{
N,
N,
N,
N};
std::array<int *, 4> tet_face_ptr{
&tet_alpha_face[0](0, 0), &tet_alpha_face[1](0, 0),
&tet_alpha_face[2](0, 0), &tet_alpha_face[3](0, 0)};
CHKERR BernsteinBezier::generateIndicesTriTet(tet_face_n.data(),
tet_face_ptr.data());
CHKERR BernsteinBezier::generateIndicesTetTet(
N, &tet_alpha_tet(0, 0));
std::vector<MatrixInt *> alpha_tet;
alpha_tet.push_back(&tet_alpha_vec);
for (int e = 0; e != 6; ++e)
alpha_tet.push_back(&tet_alpha_edge[e]);
for (
int f = 0;
f != 4; ++
f)
alpha_tet.push_back(&tet_alpha_face[f]);
alpha_tet.push_back(&tet_alpha_tet);
auto create_tet_mesh = [&](
auto &moab_ref,
Range &tets) {
std::array<double, 12> base_coords{0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
for (int nn = 0; nn < 4; nn++)
CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
const int max_level = 3;
for (int ll = 0; ll != max_level; ll++) {
CHKERR m_ref->addVerticesInTheMiddleOfEdges(
}
->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
tets);
if (1) {
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
CHKERR moab_ref.add_entities(meshset, tets);
CHKERR moab_ref.convert_entities(meshset,
true,
false,
false);
CHKERR moab_ref.delete_entities(&meshset, 1);
}
};
CHKERR create_tet_mesh(moab, tets);
CHKERR moab.get_connectivity(tets, elem_nodes,
false);
CHKERR moab.get_coords(elem_nodes, &coords(0, 0));
auto calc_lambda_on_tet = [](auto &coords) {
for (
size_t i = 0;
i != coords.size1(); ++
i) {
lambda(
i, 0) = 1 - coords(
i, 0) - coords(
i, 1) - coords(
i, 2);
}
};
auto lambda = calc_lambda_on_tet(coords);
int nn = 0;
for (auto alpha_ptr : alpha_tet) {
MatrixDouble diff_base(coords.size1(), 3 * alpha_ptr->size1());
CHKERR BernsteinBezier::baseFunctionsTet(
N, coords.size1(), alpha_ptr->size1(), &(*alpha_ptr)(0, 0),
&
lambda(0, 0), Tools::diffShapeFunMBTET.data(), &base(0, 0),
&diff_base(0, 0));
CHKERR check_property_one(coords.size1(), *alpha_ptr,
lambda, base,
BernsteinBezier::baseFunctionsTet, false);
for (
int j = 0;
j != base.size2(); ++
j) {
double def_val[] = {0};
("base_" + boost::lexical_cast<std::string>(nn) + "_" +
boost::lexical_cast<std::string>(
j))
.c_str(),
1, MB_TYPE_DOUBLE, th, MB_TAG_CREAT | MB_TAG_DENSE, def_val);
CHKERR moab.tag_set_data(th, elem_nodes, &trans_base(
j, 0));
}
++nn;
}
CHKERR moab.create_meshset(MESHSET_SET, meshset);
CHKERR moab.add_entities(meshset, tets);
CHKERR moab.write_file(
"bb.vtk",
"VTK",
"", &meshset, 1);
}
}
#define CATCH_ERRORS
Catch errors.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
UBlasMatrix< int > MatrixInt
implementation of Data Operators for Forces and Sources
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_1D_TABLE[]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Mesh refinement interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.