177 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
double *
L,
178 double *diffL,
const int dim) =
181 int nb_gauss_pts = pts.size2();
188 vert_dat.getN(base).resize(nb_gauss_pts, 6,
false);
189 vert_dat.getDiffN(base).resize(nb_gauss_pts, 18);
191 noalias(vert_dat.getDiffN(base)) =
194 auto &vert_n = vert_dat.getN(base);
195 auto &vert_diff_n = vert_dat.getDiffN(base);
197 constexpr
int prism_edge_map[9][2] = {{0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 4},
198 {2, 5}, {3, 4}, {4, 5}, {5, 3}};
200 auto edge_through_thickness = [&](
const int ee) {
208 int order = thickness_ent.getOrder();
210 if ((
unsigned int)nb_dofs != thickness_ent.getN(base).size2())
212 "nb_dofs != nb_dofs %d != %d", nb_dofs,
213 thickness_ent.getN(base).size2());
214 prism_ent.getN(base).resize(nb_gauss_pts, nb_dofs,
false);
215 prism_ent.getDiffN(base).resize(nb_gauss_pts, 3 * nb_dofs,
false);
221 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
223 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
224 double extrude = vert_n(gg, prism_edge_map[ee][0]) +
225 vert_n(gg, prism_edge_map[ee][1]);
227 double diff_extrude[3];
228 for (
auto d : {0, 1, 2})
229 diff_extrude[
d] = vert_diff_n(gg, 3 * prism_edge_map[ee][0] +
d) +
230 vert_diff_n(gg, 3 * prism_edge_map[ee][1] +
d);
232 double n0 = vert_n(gg, 0) + vert_n(gg, 1) + vert_n(gg, 2);
233 double n1 = vert_n(gg, 3) + vert_n(gg, 4) + vert_n(gg, 5);
254 double *
l = &(thickness_ent.getN(base)(ggt, 0));
255 double *diff_l_1d = &(thickness_ent.getDiffN(base)(ggt, 0));
257 double bubble = n0 * n1;
258 double diff_bubble[3];
259 for (
auto d : {0, 1, 2})
261 n0 * (vert_diff_n(gg, 3 * 3 +
d) + vert_diff_n(gg, 3 * 4 +
d) +
262 vert_diff_n(gg, 3 * 5 +
d)) +
264 n1 * (vert_diff_n(gg, 3 * 0 +
d) + vert_diff_n(gg, 3 * 1 +
d) +
265 vert_diff_n(gg, 3 * 2 +
d));
267 for (
int dd = 0;
dd != nb_dofs;
dd++) {
269 prism_ent.getN(base)(gg,
dd) = extrude * bubble *
l[
dd];
270 for (
auto d : {0, 1, 2})
271 prism_ent.getDiffN(base)(gg, 3 *
dd +
d) =
272 diff_extrude[
d] * bubble *
l[
dd] +
273 extrude * diff_bubble[
d] *
l[
dd];
275 prism_ent.getDiffN(base)(gg, 3 *
dd + 2) +=
276 extrude * bubble * 2 * diff_l_1d[
dd];
284 auto edge_on_the_triangle = [&](
const int ee) {
290 data.
dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
292 data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
294 for (
int dd = 0;
dd < nb_dofs;
dd++) {
296 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
305 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
307 double dzeta, edge_shape;
318 dksi_tri_n * edge_shape;
320 deta_tri_n * edge_shape;
329 auto trinagle_through_thickness = [&](
const int ff) {
334 data.
dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
336 data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
338 for (
int dd = 0;
dd < nb_dofs;
dd++) {
340 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
343 double dksi_tri_n[2];
344 for (
int kk = 0; kk < 2; kk++) {
349 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
351 double dzeta, edge_shape;
361 for (
auto kk : {0, 1})
363 dksi_tri_n[kk] * edge_shape;
372 auto quads_base = [&]() {
374 int quads_nodes[3 * 4];
375 int quad_order[3] = {0, 0, 0};
376 double *quad_n[3], *diff_quad_n[3];
379 SideNumber_multiIndex::nth_index<1>::type::iterator siit;
380 siit = side_table.get<1>().lower_bound(boost::make_tuple(MBQUAD, 0));
381 SideNumber_multiIndex::nth_index<1>::type::iterator hi_siit;
382 hi_siit = side_table.get<1>().upper_bound(boost::make_tuple(MBQUAD, 3));
383 if (std::distance(siit, hi_siit) != 3)
385 "Expected three quads on prisms");
389 CHKERR cTx->
mOab.get_connectivity(ent, conn, num_nodes,
true);
390 for (; siit != hi_siit; ++siit) {
391 int side = siit->get()->side_number;
392 if(side < 0 && side > 2)
394 "Side in range 0,1,2 expected");
398 CHKERR cTx->
mOab.get_connectivity(quad, conn_quad, num_nodes_quad,
true);
399 for (
int nn = 0; nn < num_nodes_quad; nn++) {
400 quads_nodes[4 * side + nn] =
401 std::distance(conn, std::find(conn, conn + 6, conn_quad[nn]));
404 quad_order[side] =
order;
413 &*data.
dataOnEntities[MBQUAD][side].getDiffN(base).data().begin();
416 diff_quad_n[side] = NULL;
419 if (quad_order[0] > 0 || quad_order[1] > 0 || quad_order[2] > 0) {
422 double *diff_vertex_n =
425 diff_vertex_n, quad_n, diff_quad_n,
426 nb_gauss_pts, base_polynomials);
431 auto prim_base = [&]() {
434 double *vertex_n = &data.
dataOnEntities[MBVERTEX][0].getN(base)(0, 0);
435 double *diff_vertex_n =
443 order, vertex_n, diff_vertex_n,
445 &data.
dataOnEntities[MBPRISM][0].getDiffN(base)(0, 0), nb_gauss_pts,
453 for (; ee <= 2; ++ee)
454 CHKERR edge_on_the_triangle(ee);
455 for (; ee <= 5; ++ee)
456 CHKERR edge_through_thickness(ee);
457 for (; ee <= 8; ++ee)
458 CHKERR edge_on_the_triangle(ee);
462 for (
int ff = 3; ff <= 4; ++ff)
463 CHKERR trinagle_through_thickness(ff);