170 {
172
175
176 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
double *
L,
177 double *diffL, const int dim) =
179
180 int nb_gauss_pts = pts.size2();
183
184
185
186 auto &vert_dat = data.dataOnEntities[MBVERTEX][0];
187 vert_dat.getN(base).resize(nb_gauss_pts, 6, false);
188 vert_dat.getDiffN(base).resize(nb_gauss_pts, 18);
189 noalias(vert_dat.getN(base)) = data.dataOnEntities[MBVERTEX][0].getN(
NOBASE);
190 noalias(vert_dat.getDiffN(base)) =
191 data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE);
192
193 auto &vert_n = vert_dat.getN(base);
194 auto &vert_diff_n = vert_dat.getDiffN(base);
195
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}};
198
199 auto edge_through_thickness = [&](const int ee) {
201
202
203
205 auto &prism_ent = data.dataOnEntities[MBEDGE][ee];
206
207 int order = thickness_ent.getOrder();
209 if ((unsigned int)nb_dofs != thickness_ent.getN(base).size2())
211 "nb_dofs != nb_dofs %d != %ld", 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);
215 if (nb_dofs == 0)
217
218
219 int gg = 0;
220 for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
221
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]);
225
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);
230
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);
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253 double *
l = &(thickness_ent.getN(base)(ggt, 0));
254 double *diff_l_1d = &(thickness_ent.getDiffN(base)(ggt, 0));
255
256 double bubble = n0 * n1;
257 double diff_bubble[3];
258 for (auto d : {0, 1, 2})
259 diff_bubble[d] =
260 n0 * (vert_diff_n(gg, 3 * 3 + d) + vert_diff_n(gg, 3 * 4 + d) +
261 vert_diff_n(gg, 3 * 5 + d)) +
262
263 n1 * (vert_diff_n(gg, 3 * 0 + d) + vert_diff_n(gg, 3 * 1 + d) +
264 vert_diff_n(gg, 3 * 2 + d));
265
266 for (
int dd = 0;
dd != nb_dofs;
dd++) {
267
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];
273
274 prism_ent.getDiffN(base)(gg, 3 *
dd + 2) +=
275 extrude * bubble * 2 * diff_l_1d[dd];
276
277 }
278 }
279 }
281 };
282
283 auto edge_on_the_triangle = [&](const int ee) {
285
286
287 int nb_dofs =
289 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
290 false);
291 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
292 3 * nb_dofs, false);
293 for (
int dd = 0;
dd < nb_dofs;
dd++) {
294 int gg = 0;
295 for (int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
298 double dksi_tri_n =
301 double deta_tri_n =
304 for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
306 double dzeta, edge_shape;
307 if (ee < 3) {
310 } else {
313 }
314 data.dataOnEntities[MBEDGE][ee].getN(base)(gg,
dd) =
315 tri_n * edge_shape;
316 data.dataOnEntities[MBEDGE][ee].getDiffN(base)(gg, 3 *
dd + 0) =
317 dksi_tri_n * edge_shape;
318 data.dataOnEntities[MBEDGE][ee].getDiffN(base)(gg, 3 *
dd + 1) =
319 deta_tri_n * edge_shape;
320 data.dataOnEntities[MBEDGE][ee].getDiffN(base)(gg, 3 *
dd + 2) =
321 tri_n * dzeta;
322 }
323 }
324 }
326 };
327
328 auto trinagle_through_thickness = [&](const int ff) {
330 int nb_dofs;
331 nb_dofs =
333 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
334 false);
335 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
336 3 * nb_dofs, false);
337 for (
int dd = 0;
dd < nb_dofs;
dd++) {
338 int gg = 0;
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++) {
344 dksi_tri_n[kk] =
347 }
348 for (int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
350 double dzeta, edge_shape;
351 if (ff == 3) {
354 } else {
357 }
358 data.dataOnEntities[MBTRI][ff].getN(base)(gg,
dd) =
359 tri_n * edge_shape;
360 for (auto kk : {0, 1})
361 data.dataOnEntities[MBTRI][ff].getDiffN(base)(gg, 3 *
dd + kk) =
362 dksi_tri_n[kk] * edge_shape;
363 data.dataOnEntities[MBTRI][ff].getDiffN(base)(gg, 3 *
dd + 2) =
364 tri_n * dzeta;
365 }
366 }
367 }
369 };
370
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");
386 int num_nodes;
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");
394 int num_nodes_quad;
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]));
401 }
402 int order = data.dataOnEntities[MBQUAD][side].getOrder();
403 quad_order[side] =
order;
404 data.dataOnEntities[MBQUAD][side].getN(base).resize(
406 data.dataOnEntities[MBQUAD][side].getDiffN(base).resize(
408 if (data.dataOnEntities[MBQUAD][side].getN(base).size2() > 0) {
409 quad_n[side] =
410 &*data.dataOnEntities[MBQUAD][side].getN(base).data().begin();
411 diff_quad_n[side] =
412 &*data.dataOnEntities[MBQUAD][side].getDiffN(base).data().begin();
413 } else {
414 quad_n[side] = NULL;
415 diff_quad_n[side] = NULL;
416 }
417 }
418 if (quad_order[0] > 0 || quad_order[1] > 0 || quad_order[2] > 0) {
419 double *vertex_n =
420 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin();
421 double *diff_vertex_n =
422 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin();
424 diff_vertex_n, quad_n, diff_quad_n,
425 nb_gauss_pts, base_polynomials);
426 }
428 };
429
430 auto prim_base = [&]() {
432 int order = data.dataOnEntities[MBPRISM][0].getOrder();
433 double *vertex_n = &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0);
434 double *diff_vertex_n =
435 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0);
436 data.dataOnEntities[MBPRISM][0].getN(base).resize(
438 data.dataOnEntities[MBPRISM][0].getDiffN(base).resize(
442 order, vertex_n, diff_vertex_n,
443 &data.dataOnEntities[MBPRISM][0].getN(base)(0, 0),
444 &data.dataOnEntities[MBPRISM][0].getDiffN(base)(0, 0), nb_gauss_pts,
445 base_polynomials);
446 }
448 };
449
450
451 int ee = 0;
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);
458
459
460
461 for (int ff = 3; ff <= 4; ++ff)
462 CHKERR trinagle_through_thickness(ff);
463
464
465
467
468
470
472}
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
PetscErrorCode H1_QuadShapeFunctions_MBPRISM(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
PetscErrorCode H1_VolumeShapeFunctions_MBPRISM(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
FTensor::Index< 'l', 3 > l
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
double zeta
Viscous hardening.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
EntitiesFieldData & dataTrianglesOnly
EntitiesFieldData & dataTroughThickness
SideNumber_multiIndex & getSideNumberTable() const
EntityHandle getEnt() const