11 using namespace MoFEM;
28 dataTrianglesOnly(data_triangles_only),
29 dataTroughThickness(data_trough_thickness),
30 gaussPtsTrianglesOnly(gauss_pts_triangles_only),
31 gaussPtsThroughThickness(gauss_pts_through_thickness), mOab(moab),
35 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
51 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
59 "Pointer to element should be given "
60 "when EntPolynomialBaseCtx is constructed "
61 "(use different constructor)");
63 int nb_gauss_pts = pts.size2();
70 "Wrong dimension of pts, should be at least 3 rows with coordinates");
77 "It is assumed that base for vertices is calculated");
82 data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 6,
false);
83 data.
dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 12,
86 (
unsigned int)nb_gauss_pts)
88 "Base functions or nodes has wrong number of integration points "
134 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
double *
L,
135 double *diffL,
const int dim) =
140 for (
unsigned int ee = 3; ee <= 5; ee++) {
142 int sense = ent_data.getSense();
143 int order = ent_data.getOrder();
145 ent_data.getN(base).resize(nb_gauss_pts_through_thickness, nb_dofs,
false);
146 ent_data.getDiffN(base).resize(nb_gauss_pts_through_thickness, nb_dofs,
151 for (
int gg = 0; gg < nb_gauss_pts_through_thickness; gg++) {
161 &ent_data.getN(base)(gg, 0),
162 &ent_data.getDiffN(base)(gg, 0), 1);
176 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
double *
L,
177 double *diffL,
const int dim) =
180 int nb_gauss_pts = pts.size2();
187 vert_dat.getN(base).resize(nb_gauss_pts, 6,
false);
188 vert_dat.getDiffN(base).resize(nb_gauss_pts, 18);
190 noalias(vert_dat.getDiffN(base)) =
193 auto &vert_n = vert_dat.getN(base);
194 auto &vert_diff_n = vert_dat.getDiffN(base);
196 constexpr
int prism_edge_map[9][2] = {{0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 4},
197 {2, 5}, {3, 4}, {4, 5}, {5, 3}};
199 auto edge_through_thickness = [&](
const int ee) {
207 int order = thickness_ent.getOrder();
209 if ((
unsigned int)nb_dofs != thickness_ent.getN(base).size2())
211 "nb_dofs != nb_dofs %d != %d", nb_dofs,
212 thickness_ent.getN(base).size2());
213 prism_ent.getN(base).resize(nb_gauss_pts, nb_dofs,
false);
214 prism_ent.getDiffN(base).resize(nb_gauss_pts, 3 * nb_dofs,
false);
220 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
222 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
223 double extrude = vert_n(gg, prism_edge_map[ee][0]) +
224 vert_n(gg, prism_edge_map[ee][1]);
226 double diff_extrude[3];
227 for (
auto d : {0, 1, 2})
228 diff_extrude[
d] = vert_diff_n(gg, 3 * prism_edge_map[ee][0] +
d) +
229 vert_diff_n(gg, 3 * prism_edge_map[ee][1] +
d);
231 double n0 = vert_n(gg, 0) + vert_n(gg, 1) + vert_n(gg, 2);
232 double n1 = vert_n(gg, 3) + vert_n(gg, 4) + vert_n(gg, 5);
253 double *
l = &(thickness_ent.getN(base)(ggt, 0));
254 double *diff_l_1d = &(thickness_ent.getDiffN(base)(ggt, 0));
256 double bubble = n0 * n1;
257 double diff_bubble[3];
258 for (
auto d : {0, 1, 2})
260 n0 * (vert_diff_n(gg, 3 * 3 +
d) + vert_diff_n(gg, 3 * 4 +
d) +
261 vert_diff_n(gg, 3 * 5 +
d)) +
263 n1 * (vert_diff_n(gg, 3 * 0 +
d) + vert_diff_n(gg, 3 * 1 +
d) +
264 vert_diff_n(gg, 3 * 2 +
d));
266 for (
int dd = 0;
dd != nb_dofs;
dd++) {
268 prism_ent.getN(base)(gg,
dd) = extrude * bubble *
l[
dd];
269 for (
auto d : {0, 1, 2})
270 prism_ent.getDiffN(base)(gg, 3 *
dd +
d) =
271 diff_extrude[
d] * bubble *
l[
dd] +
272 extrude * diff_bubble[
d] *
l[
dd];
274 prism_ent.getDiffN(base)(gg, 3 *
dd + 2) +=
275 extrude * bubble * 2 * diff_l_1d[
dd];
283 auto edge_on_the_triangle = [&](
const int ee) {
289 data.
dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
291 data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
293 for (
int dd = 0;
dd < nb_dofs;
dd++) {
295 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
304 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
306 double dzeta, edge_shape;
317 dksi_tri_n * edge_shape;
319 deta_tri_n * edge_shape;
328 auto trinagle_through_thickness = [&](
const int ff) {
333 data.
dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
335 data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
337 for (
int dd = 0;
dd < nb_dofs;
dd++) {
339 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
342 double dksi_tri_n[2];
343 for (
int kk = 0; kk < 2; kk++) {
348 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
350 double dzeta, edge_shape;
360 for (
auto kk : {0, 1})
362 dksi_tri_n[kk] * edge_shape;
371 auto quads_base = [&]() {
373 int quads_nodes[3 * 4];
374 int quad_order[3] = {0, 0, 0};
375 double *quad_n[3], *diff_quad_n[3];
378 SideNumber_multiIndex::nth_index<1>::type::iterator siit;
379 siit = side_table.get<1>().lower_bound(boost::make_tuple(MBQUAD, 0));
380 SideNumber_multiIndex::nth_index<1>::type::iterator hi_siit;
381 hi_siit = side_table.get<1>().upper_bound(boost::make_tuple(MBQUAD, 3));
382 if (std::distance(siit, hi_siit) != 3)
384 "Expected three quads on prisms");
388 CHKERR cTx->
mOab.get_connectivity(ent, conn, num_nodes,
true);
389 for (; siit != hi_siit; ++siit) {
390 int side = siit->get()->side_number;
391 if(side < 0 && side > 2)
393 "Side in range 0,1,2 expected");
397 CHKERR cTx->
mOab.get_connectivity(quad, conn_quad, num_nodes_quad,
true);
398 for (
int nn = 0; nn < num_nodes_quad; nn++) {
399 quads_nodes[4 * side + nn] =
400 std::distance(conn, std::find(conn, conn + 6, conn_quad[nn]));
403 quad_order[side] =
order;
412 &*data.
dataOnEntities[MBQUAD][side].getDiffN(base).data().begin();
415 diff_quad_n[side] = NULL;
418 if (quad_order[0] > 0 || quad_order[1] > 0 || quad_order[2] > 0) {
421 double *diff_vertex_n =
424 diff_vertex_n, quad_n, diff_quad_n,
425 nb_gauss_pts, base_polynomials);
430 auto prim_base = [&]() {
433 double *vertex_n = &data.
dataOnEntities[MBVERTEX][0].getN(base)(0, 0);
434 double *diff_vertex_n =
442 order, vertex_n, diff_vertex_n,
444 &data.
dataOnEntities[MBPRISM][0].getDiffN(base)(0, 0), nb_gauss_pts,
452 for (; ee <= 2; ++ee)
453 CHKERR edge_on_the_triangle(ee);
454 for (; ee <= 5; ++ee)
455 CHKERR edge_through_thickness(ee);
456 for (; ee <= 8; ++ee)
457 CHKERR edge_on_the_triangle(ee);
461 for (
int ff = 3; ff <= 4; ++ff)
462 CHKERR trinagle_through_thickness(ff);