8static char help[] =
"testing interface inserting algorithm\n\n";
12 for (
unsigned int ii = 0; ii <
m.size1(); ii++) {
13 for (
unsigned int jj = 0; jj <
m.size2(); jj++) {
20int main(
int argc,
char *argv[]) {
26 moab::Core mb_instance;
27 moab::Interface &moab = mb_instance;
35 INTEGRATEDJACOBIPOLYNOMIAL,
37 H1TET_BERNSTEIN_BEZIER,
44 H1TRI_BERNSTEIN_BEZIER,
52 H1EDGE_BERNSTEIN_BEZIER,
60 const char *list[] = {
"legendre",
65 "h1tet_bernstein_bezier",
72 "h1tri_bernstein_bezier",
80 "h1edge_bernstein_bezier",
81 "hcurledge_ainsworth",
82 "hcurledge_demkowicz",
87 PetscInt choice_value = LEGENDREPOLYNOMIAL;
91 if (flg != PETSC_TRUE) {
100 boost::shared_ptr<MatrixDouble> base_ptr(
new MatrixDouble);
101 boost::shared_ptr<MatrixDouble> diff_base_ptr(
new MatrixDouble);
103 const double eps = 1e-3;
105 if (choice_value == LEGENDREPOLYNOMIAL) {
109 4, &diff_s, 1, base_ptr, diff_base_ptr)));
112 std::cout <<
"LegendrePolynomial\n";
113 std::cout << pts_1d << std::endl;
114 std::cout << *base_ptr << std::endl;
115 std::cout << *diff_base_ptr << std::endl;
118 std::cout << sum << std::endl;
119 std::cout << diff_sum << std::endl;
120 if (fabs(2.04688 - sum) >
eps) {
123 if (fabs(2.25 - diff_sum) >
eps) {
128 pts_1d.resize(1, 11,
false);
129 for (
int ii = 0; ii != 11; ii++) {
130 pts_1d(0, ii) = 2 * ((
double)ii / 10) - 1;
133 boost::shared_ptr<MatrixDouble> kernel_base_ptr(
new MatrixDouble);
134 boost::shared_ptr<MatrixDouble> diff_kernel_base_ptr(
new MatrixDouble);
136 if (choice_value == LOBATTOPOLYNOMIAL) {
140 7 + 2, &diff_s, 1, base_ptr, diff_base_ptr)));
145 7, &diff_s, 1, kernel_base_ptr, diff_kernel_base_ptr)));
147 for (
int ii = 0; ii != 11; ii++) {
148 double s = pts_1d(0, ii);
149 std::cerr <<
"lobatto_plot " << s <<
" " << (*base_ptr)(ii, 1) <<
" "
150 << (*diff_base_ptr)(ii, 1) <<
" " << (*kernel_base_ptr)(ii, 1)
151 <<
" " << (*diff_kernel_base_ptr)(ii, 1) <<
" "
152 << (*kernel_base_ptr)(ii, 1) * (1 - s * s) <<
" "
153 << (*kernel_base_ptr)(ii, 1) * (-2 * s) +
154 (*diff_kernel_base_ptr)(ii, 1) * (1 - s * s)
157 std::cout <<
"LobattoPolynomial\n";
158 std::cout << pts_1d << std::endl;
159 std::cout << *base_ptr << std::endl;
160 std::cout << *diff_base_ptr << std::endl;
164 std::cout << sum << std::endl;
165 std::cout << diff_sum << std::endl;
166 if (fabs(9.39466 - sum) >
eps) {
169 if (fabs(14.0774 - diff_sum) >
eps) {
175 double diff_sum =
sum_matrix(*diff_kernel_base_ptr);
176 std::cout << sum << std::endl;
177 std::cout << diff_sum << std::endl;
178 if (fabs(-13.9906 * 4 - sum) >
eps) {
181 if (fabs(-101.678 * 4 - diff_sum) >
eps) {
187 if (choice_value == JACOBIPOLYNOMIAL) {
191 pts_1d.resize(1,
n,
false);
193 for (
int ii = 0; ii !=
n; ii++) {
194 pts_1d(0, ii) = (
double)ii / (
n - 1);
199 diff_base_ptr->clear();
206 5, &diff_x, &diff_t, 1, 0, base_ptr, diff_base_ptr)));
208 for (
int ii = 0; ii !=
n; ii++) {
209 double s = pts_1d(0, ii);
210 std::cerr <<
"jacobi_plot " << s <<
" " << (*base_ptr)(ii, 4) <<
" "
211 << (*diff_base_ptr)(ii, 4) << std::endl;
213 for (
int ii = 0; ii !=
n - 1; ii++) {
214 double s = (pts_1d(0, ii) + pts_1d(0, ii + 1)) / 2;
215 std::cerr <<
"jacobi_diff_plot_fd " << s <<
" "
216 << ((*base_ptr)(ii + 1, 4) - (*base_ptr)(ii, 4)) /
220 std::cout <<
"JacobiPolynomial\n";
221 std::cout << pts_1d << std::endl;
222 std::cout << *base_ptr << std::endl;
223 std::cout << *diff_base_ptr << std::endl;
227 std::cout << sum << std::endl;
228 std::cout << diff_sum << std::endl;
229 if (fabs(23.2164 - sum) >
eps) {
232 if (fabs(167.995 - diff_sum) >
eps) {
238 if (choice_value == INTEGRATEDJACOBIPOLYNOMIAL) {
242 pts_1d.resize(1,
n,
false);
244 for (
int ii = 0; ii !=
n; ii++) {
245 pts_1d(0, ii) = (
double)ii / 20.;
246 pts_1d_t(0, ii) = 1 - pts_1d(0, ii);
250 diff_base_ptr->clear();
257 6, &diff_x, &diff_t, 1, 1, base_ptr, diff_base_ptr)));
259 for (
int ii = 0; ii !=
n; ii++) {
260 double s = pts_1d(0, ii);
261 std::cerr <<
"integrated_jacobi_plot " << s <<
" " << (
double)ii / 20.
262 <<
" " << (*base_ptr)(ii, 1) <<
" " << (*base_ptr)(ii, 2)
263 <<
" " << (*base_ptr)(ii, 3) <<
" " << (*base_ptr)(ii, 4)
264 <<
" " << (*base_ptr)(ii, 5) << endl;
268 std::cout <<
"IntegratedJacobiPolynomial\n";
269 std::cout << pts_1d << std::endl;
270 std::cout << *base_ptr << std::endl;
271 std::cout << *diff_base_ptr << std::endl;
275 std::cout << sum << std::endl;
276 std::cout << diff_sum << std::endl;
287 for (
int type = MBVERTEX; type != MBMAXTYPE; type++) {
296 for (
int ee = 0; ee < 6; ee++) {
301 for (
int ff = 0; ff < 4; ff++) {
307 tet_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(4,
false);
308 std::fill(tet_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
314 pts_tet.resize(3, nb_gauss_pts,
false);
315 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[1], 4,
317 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[2], 4,
319 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[3], 4,
326 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[tet_rule]->points, 1,
330 const int cannonical_tet_face[4][3] = {
331 {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
332 for (
int ff = 0; ff < 4; ff++) {
333 for (
int nn = 0; nn < 3; nn++) {
334 tet_data.
facesNodes(ff, nn) = cannonical_tet_face[ff][nn];
338 if (choice_value == H1TET_AINSWORTH) {
344 &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
356 "Different pointers");
359 double sum = 0, diff_sum = 0;
360 std::cout <<
"Edges\n";
361 for (
int ee = 0; ee < 6; ee++) {
373 std::cout <<
"Faces\n";
374 for (
int ff = 0; ff < 4; ff++) {
386 std::cout <<
"Tets\n";
397 std::cout <<
"sum " << sum << std::endl;
398 std::cout <<
"diff_sum " << diff_sum << std::endl;
399 if (fabs(1.3509 - sum) >
eps) {
402 if (fabs(0.233313 - diff_sum) >
eps) {
407 if (choice_value == H1TET_BERNSTEIN_BEZIER) {
413 &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
416 tet_data,
"TEST_FIELD",
H1,
424 double sum = 0, diff_sum = 0;
425 std::cout <<
"Vertices\n";
426 std::cout << tet_data.
dataOnEntities[MBVERTEX][0].getN(
"TEST_FIELD")
428 std::cout << tet_data.
dataOnEntities[MBVERTEX][0].getDiffN(
"TEST_FIELD")
434 std::cout <<
"Edges\n";
435 for (
int ee = 0; ee < 6; ee++) {
436 std::cout << tet_data.
dataOnEntities[MBEDGE][ee].getN(
"TEST_FIELD")
438 std::cout << tet_data.
dataOnEntities[MBEDGE][ee].getDiffN(
"TEST_FIELD")
445 std::cout <<
"Faces\n";
446 for (
int ff = 0; ff < 4; ff++) {
447 std::cout << tet_data.
dataOnEntities[MBTRI][ff].getN(
"TEST_FIELD")
449 std::cout << tet_data.
dataOnEntities[MBTRI][ff].getDiffN(
"TEST_FIELD")
456 std::cout <<
"Tets\n";
459 std::cout << tet_data.
dataOnEntities[MBTET][0].getDiffN(
"TEST_FIELD")
464 std::cout <<
"sum " << sum << std::endl;
465 std::cout <<
"diff_sum " << diff_sum << std::endl;
466 if (fabs(4.38395 - sum) >
eps) {
469 if (fabs(diff_sum) >
eps) {
474 if (choice_value == HDIVTET_AINSWORTH) {
479 double sum = 0, diff_sum = 0;
480 std::cout <<
"Faces\n";
481 for (
int ff = 0; ff < 4; ff++) {
493 std::cout <<
"Tets\n";
504 std::cout <<
"sum " << 1e8 * sum << std::endl;
505 std::cout <<
"diff_sum " << 1e8 * diff_sum << std::endl;
506 if (fabs(0.188636 - sum) >
eps) {
509 if (fabs(32.9562 - diff_sum) >
eps) {
514 if (choice_value == HDIVTET_DEMKOWICZ) {
519 double sum = 0, diff_sum = 0;
520 std::cout <<
"Faces\n";
521 for (
int ff = 0; ff < 4; ff++) {
533 std::cout <<
"Tets\n";
543 std::cout <<
"sum " << sum << std::endl;
544 std::cout <<
"diff_sum " << diff_sum << std::endl;
545 const double expected_result = -2.70651;
546 const double expected_diff_result = 289.421;
547 if (fabs((expected_result - sum) / expected_result) >
eps) {
550 if (fabs((expected_diff_result - diff_sum) / expected_diff_result) >
556 if (choice_value == HCURLTET_AINSWORTH) {
561 double sum = 0, diff_sum = 0;
562 std::cout <<
"Edges\n";
563 for (
int ee = 0; ee < 6; ee++) {
575 std::cout <<
"Faces\n";
576 for (
int ff = 0; ff < 4; ff++) {
588 std::cout <<
"Tets\n";
599 std::cout <<
"sum " << sum << std::endl;
600 std::cout <<
"diff_sum " << diff_sum << std::endl;
601 if (fabs(-1.7798 - sum) >
eps) {
604 if (fabs(-67.1793 - diff_sum) >
eps) {
609 if (choice_value == HCURLTET_DEMKOWICZ) {
613 double sum = 0, diff_sum = 0;
614 std::cout <<
"Edges\n";
615 for (
int ee = 0; ee < 6; ee++) {
627 std::cout <<
"Faces\n";
628 for (
int ff = 0; ff < 4; ff++) {
640 std::cout <<
"Tets\n";
650 std::cout <<
"sum " << sum << std::endl;
651 std::cout <<
"diff_sum " << diff_sum << std::endl;
652 if (fabs(7.35513 - sum) >
eps) {
655 if (fabs(62.4549 - diff_sum) >
eps) {
660 if (choice_value == L2TET) {
665 double sum = 0, diff_sum = 0;
666 std::cout <<
"Tets\n";
677 std::cout <<
"sum " << sum << std::endl;
678 std::cout <<
"diff_sum " << diff_sum << std::endl;
679 if (fabs(3.60352 - sum) >
eps) {
682 if (fabs(-36.9994 - diff_sum) >
eps) {
688 for (
int type = MBVERTEX; type != MBMAXTYPE; type++) {
697 for (
int ee = 0; ee < 3; ee++) {
703 tri_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(3,
false);
704 std::fill(tri_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
710 pts_tri.resize(2, nb_gauss_pts,
false);
711 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[tri_rule]->points[1], 3,
713 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[tri_rule]->points[2], 3,
720 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[tri_rule]->points, 1,
728 if (choice_value == H1TRI_AINSWORTH) {
737 "Different pointers");
739 double sum = 0, diff_sum = 0;
740 std::cout <<
"Edges\n";
741 for (
int ee = 0; ee < 3; ee++) {
753 std::cout <<
"Face\n";
764 std::cout <<
"sum " << sum << std::endl;
765 std::cout <<
"diff_sum " << diff_sum << std::endl;
766 if (fabs(0.805556 - sum) >
eps)
769 if (fabs(0.0833333 - diff_sum) >
eps)
773 if (choice_value == H1TRI_BERNSTEIN_BEZIER) {
776 tri_data,
"TET_FIELD",
H1,
783 "Pointers should be diffrent");
785 double sum = 0, diff_sum = 0;
786 std::cout <<
"Vertex\n";
787 std::cout << tri_data.
dataOnEntities[MBVERTEX][0].getN(
"TET_FIELD")
789 std::cout << tri_data.
dataOnEntities[MBVERTEX][0].getDiffN(
"TET_FIELD")
794 std::cout <<
"Edges\n";
795 for (
int ee = 0; ee < 3; ee++) {
796 std::cout << tri_data.
dataOnEntities[MBEDGE][ee].getN(
"TET_FIELD")
798 std::cout << tri_data.
dataOnEntities[MBEDGE][ee].getDiffN(
"TET_FIELD")
805 std::cout <<
"Face\n";
808 std::cout << tri_data.
dataOnEntities[MBTRI][0].getDiffN(
"TET_FIELD")
813 std::cout <<
"sum " << sum << std::endl;
814 std::cout <<
"diff_sum " << diff_sum << std::endl;
815 if (std::abs(3.01389 - sum) >
eps)
818 if (std::abs(diff_sum) >
eps)
820 "wrong result %3.4e != $3.4e", 0, diff_sum);
823 if (choice_value == HDIVTRI_AINSWORTH) {
832 "Different pointers");
835 std::cout <<
"Face\n";
841 std::cout <<
"sum " << sum << std::endl;
842 if (fabs(1.93056 - sum) >
eps)
846 if (choice_value == HDIVTRI_DEMKOWICZ) {
856 "Different pointers");
859 std::cout <<
"Face\n";
864 std::cout <<
"sum " << sum << std::endl;
865 const double expected_result = 28.25;
866 if (fabs((expected_result - sum) / expected_result) >
eps) {
871 if (choice_value == HCURLTRI_AINSWORTH) {
881 "Different pointers");
884 std::cout <<
"Edges\n";
885 for (
int ee = 0; ee < 3; ee++) {
892 std::cout <<
"Face\n";
898 std::cout <<
"sum " << sum << std::endl;
899 if (fabs(0.333333 - sum) >
eps) {
904 if (choice_value == HCURLTRI_DEMKOWICZ) {
913 "Different pointers");
916 std::cout <<
"Edges\n";
917 for (
int ee = 0; ee < 3; ee++) {
924 std::cout <<
"Face\n";
929 std::cout <<
"sum " << sum << std::endl;
930 if (fabs(1.04591 - sum) >
eps) {
935 if (choice_value == L2TRI) {
945 "Different pointers");
948 std::cout <<
"Face\n";
954 std::cout <<
"sum " << sum << std::endl;
955 if (fabs(0.671875 - sum) >
eps) {
961 for (
int type = MBVERTEX; type != MBMAXTYPE; type++) {
972 edge_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(2,
false);
973 std::fill(edge_data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
979 pts_edge.resize(1, nb_gauss_pts,
false);
980 cblas_dcopy(nb_gauss_pts, &
QUAD_1D_TABLE[edge_rule]->points[1], 2,
987 cblas_dcopy(2 * nb_gauss_pts,
QUAD_1D_TABLE[edge_rule]->points, 1,
991 if (choice_value == H1EDGE_AINSWORTH) {
1000 "Different pointers");
1002 double sum = 0, diff_sum = 0;
1003 std::cout <<
"Edge\n";
1014 std::cout <<
"sum " << sum << std::endl;
1015 std::cout <<
"diff sum " << diff_sum << std::endl;
1016 if (std::abs(0.506122 - sum) >
eps)
1019 if (std::abs(2.85714 - diff_sum) >
eps)
1023 if (choice_value == H1EDGE_BERNSTEIN_BEZIER) {
1026 edge_data,
"TET_FIELD",
H1,
1033 "Should be diffrent pointers");
1035 double sum = 0, diff_sum = 0;
1036 std::cout <<
"Vertex\n";
1037 std::cout << edge_data.
dataOnEntities[MBVERTEX][0].getN(
"TET_FIELD")
1039 std::cout << edge_data.
dataOnEntities[MBVERTEX][0].getDiffN(
"TET_FIELD")
1041 std::cout <<
"Edge\n";
1042 std::cout << edge_data.
dataOnEntities[MBEDGE][0].getN(
"TET_FIELD")
1044 std::cout << edge_data.
dataOnEntities[MBEDGE][0].getDiffN(
"TET_FIELD")
1053 std::cout <<
"sum " << sum << std::endl;
1054 std::cout <<
"diff sum " << diff_sum << std::endl;
1055 if (std::abs(4 - sum) >
eps)
1058 if (std::abs(diff_sum) >
eps)
1062 if (choice_value == HCURLEDGE_AINSWORTH) {
1072 "Different pointers");
1080 std::cout <<
"sum " << sum << std::endl;
1081 if (fabs(-4 - sum) >
eps) {
1086 if (choice_value == HCURLEDGE_DEMKOWICZ) {
1095 "Different pointers");
1103 std::cout <<
"sum " << sum << std::endl;
1104 if (fabs(4 - sum) >
eps) {
1110 for (
int type = MBVERTEX; type != MBMAXTYPE; type++) {
1116 for (
int ee = 0; ee < 4; ee++) {
1128 nb_gauss_pts = pts_quad.size2();
1133 for (
int i = 0;
i < nb_gauss_pts; ++
i) {
1161 if (choice_value == H1QUAD) {
1170 "Different pointers");
1172 double sum = 0, diff_sum = 0;
1173 for (
int ee = 0; ee < 4; ee++) {
1184 std::cout << sum <<
" " << diff_sum << endl;
1186 if (std::abs(3.58249 - sum) >
eps)
1190 if (std::abs(-0.134694 - diff_sum) >
eps)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_1D_TABLE[]
static QUAD *const QUAD_3D_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.
Calculate base functions on tetrahedral.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to pass element data to calculate base functions on tet,triangle,edge.
data structure for finite element entity
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Legendre base functions.
Calculating Legendre base functions.
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Kernel Lobatto base functions.
Calculating Lobatto base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Legendre base functions.
Calculating Legendre base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Lobatto base functions.
Calculating Lobatto base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on triangle.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on tetrahedral.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on triangle.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
static double sum_matrix(MatrixDouble &m)