20 #include "libmesh/cell_pyramid13.h" 21 #include "libmesh/edge_edge3.h" 22 #include "libmesh/face_tri6.h" 23 #include "libmesh/face_quad8.h" 24 #include "libmesh/enum_io_package.h" 25 #include "libmesh/enum_order.h" 41 {0, 1, 4, 5, 10, 9, 99, 99},
42 {1, 2, 4, 6, 11, 10, 99, 99},
43 {2, 3, 4, 7, 12, 11, 99, 99},
44 {3, 0, 4, 8, 9, 12, 99, 99},
45 {0, 3, 2, 1, 8, 7, 6, 5}
89 const unsigned int s)
const 91 libmesh_assert_less (s,
n_sides());
100 libmesh_assert_less(s,
n_sides());
101 auto trim = (s == 4) ? 0 : 2;
105 std::vector<unsigned>
108 libmesh_assert_less(e,
n_edges());
113 const unsigned int e)
const 115 libmesh_assert_less (e,
n_edges());
140 unsigned int side_node)
const 142 libmesh_assert_less (side, this->
n_sides());
156 unsigned int edge_node)
const 158 libmesh_assert_less(edge, this->
n_edges());
168 libmesh_assert_less (i, this->
n_sides());
170 std::unique_ptr<Elem> face;
179 face = std::make_unique<Tri6>();
184 face = std::make_unique<Quad8>();
188 libmesh_error_msg(
"Invalid side i = " << i);
192 for (
auto n : face->node_index_range())
195 face->set_interior_parent(
this);
196 face->inherit_data_from(*
this);
204 const unsigned int i)
206 libmesh_assert_less (i, this->
n_sides());
215 if (!side.get() || side->type() !=
TRI6)
224 if (!side.get() || side->type() !=
QUAD8)
232 libmesh_error_msg(
"Invalid side i = " << i);
235 side->inherit_data_from(*
this);
238 for (
auto n : side->node_index_range())
246 return this->simple_build_edge_ptr<Edge3,Pyramid13>(i);
253 this->simple_build_edge_ptr<Pyramid13>(edge, i,
EDGE3);
260 std::vector<dof_id_type> & )
const 271 libmesh_not_implemented();
277 libmesh_not_implemented();
281 libmesh_error_msg(
"Unsupported IO package " << iop);
302 libmesh_error_msg(
"Invalid node n = " << n);
308 const unsigned int v)
const 310 libmesh_assert_greater_equal (n, this->
n_vertices());
311 libmesh_assert_less (n, this->
n_nodes());
324 libmesh_assert_less (v, 2);
329 unsigned short node_list[8][2] =
341 return node_list[n-5][v];
345 libmesh_error_msg(
"Invalid n = " << n);
369 x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
370 -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
371 x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
372 x0/4 - x1/4 + x2/4 - x3/4,
373 -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
374 x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
375 -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
376 x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
377 x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
378 -x0 - x1 - x2 - x3 + 2*x5 + 2*x7,
379 x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
380 -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
381 x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7
389 x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
390 -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
391 x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
392 x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
393 -x0 - x1 - x2 - x3 + 2*x6 + 2*x8,
394 x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
395 x0/4 - x1/4 + x2/4 - x3/4,
396 -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
397 x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
398 -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
399 x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
400 -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
401 x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
408 x0/4 + x1/4 + x10 + x11 + x12 + x2/4 + x3/4 - x4 - x5 - x6 - x7 - x8 + x9,
409 -3*x0/4 - 3*x1/4 - 5*x10 - 5*x11 - 5*x12 - 3*x2/4 - 3*x3/4 + 7*x4 + 4*x5 + 4*x6 + 4*x7 + 4*x8 - 5*x9,
410 3*x0/4 + 3*x1/4 + 9*x10 + 9*x11 + 9*x12 + 3*x2/4 + 3*x3/4 - 15*x4 - 6*x5 - 6*x6 - 6*x7 - 6*x8 + 9*x9,
411 -x0/4 - x1/4 - 7*x10 - 7*x11 - 7*x12 - x2/4 - x3/4 + 13*x4 + 4*x5 + 4*x6 + 4*x7 + 4*x8 - 7*x9,
412 2*x10 + 2*x11 + 2*x12 - 4*x4 - x5 - x6 - x7 - x8 + 2*x9,
413 x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
414 -3*x0/4 - 3*x1/4 + 3*x10 - 3*x11 - 3*x12 + 3*x2/4 + 3*x3/4 - 3*x5/2 + 3*x7/2 + 3*x9,
415 3*x0/4 + 3*x1/4 - 3*x10 + 3*x11 + 3*x12 - 3*x2/4 - 3*x3/4 + 3*x5/2 - 3*x7/2 - 3*x9,
416 -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
417 x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
418 -3*x0/4 + 3*x1/4 - 3*x10 - 3*x11 + 3*x12 + 3*x2/4 - 3*x3/4 + 3*x6/2 - 3*x8/2 + 3*x9,
419 3*x0/4 - 3*x1/4 + 3*x10 + 3*x11 - 3*x12 - 3*x2/4 + 3*x3/4 - 3*x6/2 + 3*x8/2 - 3*x9,
420 -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
421 -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
422 x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
423 -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
424 x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
425 -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
426 x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
430 static const int dx_dxi_exponents[14][3] =
449 static const int dx_deta_exponents[14][3] =
468 static const int dx_dzeta_exponents[19][3] =
497 a1 = -7.1805574131988893873307823958101e-01_R,
498 a2 = -5.0580870785392503961340276902425e-01_R,
499 a3 = -2.2850430565396735359574351631314e-01_R,
501 b1 = 7.2994024073149732155837979012003e-02_R,
502 b2 = 3.4700376603835188472176354340395e-01_R,
503 b3 = 7.0500220988849838312239847758405e-01_R,
506 w1 = 4.8498876871878584357513834016440e-02_R,
507 w2 = 4.5137737425884574692441981593901e-02_R,
508 w3 = 9.2440441384508327195915094925393e-03_R,
509 w4 = 7.7598202995005734972022134426305e-02_R,
510 w5 = 7.2220379881415319507907170550242e-02_R,
511 w6 = 1.4790470621521332351346415188063e-02_R,
512 w7 = 1.2415712479200917595523541508209e-01_R,
513 w8 = 1.1555260781026451121265147288039e-01_R,
514 w9 = 2.3664752994434131762154264300901e-02_R;
517 static const Real xi[N][3] =
548 static const Real eta[N][3] =
579 static const Real zeta[N][5] =
581 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
582 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
583 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
584 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
585 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
586 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
587 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
588 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
589 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
590 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
591 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
592 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
593 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
594 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
595 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
596 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
597 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
598 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
599 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
600 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
601 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
602 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
603 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
604 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3},
605 { 1., b1, b1*b1, b1*b1*b1, b1*b1*b1*b1},
606 { 1., b2, b2*b2, b2*b2*b2, b2*b2*b2*b2},
607 { 1., b3, b3*b3, b3*b3*b3, b3*b3*b3*b3}
610 static const Real w[N] = {w1, w2, w3, w4, w5, w6,
611 w1, w2, w3, w4, w5, w6,
612 w7, w8, w9, w4, w5, w6,
613 w1, w2, w3, w4, w5, w6,
617 for (
int q=0; q<N; ++q)
621 den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
622 den3 = den2*(1. - zeta[q][1]);
625 Point dx_dxi_q, dx_deta_q;
626 for (
int c=0; c<14; ++c)
629 xi[q][dx_dxi_exponents[c][0]]*
630 eta[q][dx_dxi_exponents[c][1]]*
631 zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];
634 xi[q][dx_deta_exponents[c][0]]*
635 eta[q][dx_deta_exponents[c][1]]*
636 zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
641 for (
int c=0; c<19; ++c)
644 xi[q][dx_dzeta_exponents[c][0]]*
645 eta[q][dx_dzeta_exponents[c][1]]*
646 zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
664 libmesh_assert_less (perm_num, 4);
666 for (
unsigned int i = 0; i != perm_num; ++i)
695 libmesh_assert_less (s, 5);
ElemType
Defines an enum for geometric element types.
void swap2boundaryedges(unsigned short e1, unsigned short e2, BoundaryInfo *boundary_info) const
Swaps two edges in boundary_info, if it is non-null.
static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]
This maps the node of the edge to element node numbers.
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
This maps the node of the side to element node numbers.
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
virtual bool has_affine_map() const override
Order
defines an enum for polynomial orders.
Node ** _nodes
Pointers to the nodes we are connected to.
virtual unsigned int n_nodes() const override
virtual Real volume() const override
Specialization for computing the volume of a Pyramid13.
static const int num_edges
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
virtual void connectivity(const unsigned int sc, const IOPackage iop, std::vector< dof_id_type > &conn) const override
virtual std::vector< unsigned int > nodes_on_edge(const unsigned int e) const override
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const override
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
static const int nodes_per_edge
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i) override
Builds a EDGE3 coincident with edge i.
void swap2boundarysides(unsigned short s1, unsigned short s2, BoundaryInfo *boundary_info) const
Swaps two sides in boundary_info, if it is non-null.
virtual unsigned int n_sides() const override
static const int nodes_per_side
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int n_vertices() const override
virtual bool is_edge(const unsigned int i) const override
virtual bool is_vertex(const unsigned int i) const override
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
void swap4nodes(unsigned int n1, unsigned int n2, unsigned int n3, unsigned int n4)
Swaps four node_ptrs, "rotating" them.
ElemMappingType mapping_type() const
void swap2nodes(unsigned int n1, unsigned int n2)
Swaps two node_ptrs.
virtual unsigned int n_sub_elem() const override
FIXME: we don't yet have a refinement pattern for pyramids...
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
static const int num_nodes
Geometric constants for Pyramid13.
virtual void flip(BoundaryInfo *) override final
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
static const int num_sides
Geometric constants for all Pyramids.
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i) override
Builds a QUAD8 or TRI6 coincident with face i.
void swap2neighbors(unsigned int n1, unsigned int n2)
Swaps two neighbor_ptrs.
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const override
virtual bool is_face(const unsigned int i) const override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
virtual unsigned int n_edges() const override
virtual void permute(unsigned int perm_num) override final
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
ElemType side_type(const unsigned int s) const override final
virtual Real volume() const
void swap4neighbors(unsigned int n1, unsigned int n2, unsigned int n3, unsigned int n4)
Swaps four neighbor_ptrs, "rotating" them.
virtual Order default_order() const override
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const override
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Point & point(const unsigned int i) const
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const override