24 MatrixDouble &gauss_pts_through_thickness, moab::Interface &moab,
29 dataTrianglesOnly(data_triangles_only),
30 dataTroughThickness(data_trough_thickness),
31 gaussPtsTrianglesOnly(gauss_pts_triangles_only),
32 gaussPtsThroughThickness(gauss_pts_through_thickness), mOab(moab),
36 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
52 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
60 "Pointer to element should be given "
61 "when EntPolynomialBaseCtx is constructed "
62 "(use different constructor)");
64 int nb_gauss_pts = pts.size2();
71 "Wrong dimension of pts, should be at least 3 rows with coordinates");
78 "It is assumed that base for vertices is calculated");
83 data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 6,
false);
84 data.
dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 12,
87 (
unsigned int)nb_gauss_pts)
89 "Base functions or nodes has wrong number of integration points "
135 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
double *
L,
136 double *diffL,
const int dim) =
141 for (
unsigned int ee = 3; ee <= 5; ee++) {
143 int sense = ent_data.getSense();
144 int order = ent_data.getOrder();
146 ent_data.getN(base).resize(nb_gauss_pts_through_thickness, nb_dofs,
false);
147 ent_data.getDiffN(base).resize(nb_gauss_pts_through_thickness, nb_dofs,
152 for (
int gg = 0; gg < nb_gauss_pts_through_thickness; gg++) {
162 &ent_data.getN(base)(gg, 0),
163 &ent_data.getDiffN(base)(gg, 0), 1);
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);
FieldApproximationBase
approximation base
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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)
Numer 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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
double zeta
Viscous hardening.
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
const FieldApproximationBase bAse
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Class used to pass element data to calculate base functions on fat prism.
const NumeredEntFiniteElement * fePtr
FatPrismPolynomialBaseCtx(EntitiesFieldData &data, EntitiesFieldData &data_triangles_only, EntitiesFieldData &data_trough_thickness, MatrixDouble &gauss_pts_triangles_only, MatrixDouble &gauss_pts_through_thickness, moab::Interface &moab, const NumeredEntFiniteElement *fe_ptr, const FieldSpace space, const FieldApproximationBase base, const FieldApproximationBase copy_node_base=LASTBASE)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MatrixDouble & gaussPtsThroughThickness
EntitiesFieldData & dataTrianglesOnly
EntitiesFieldData & dataTroughThickness
~FatPrismPolynomialBaseCtx()
MatrixDouble & gaussPtsTrianglesOnly
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
FatPrismPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
MoFEMErrorCode getValueH1ThroughThickness()
MoFEMErrorCode getValueH1TrianglesOnly()
~FatPrismPolynomialBase()
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Class used to pass element data to calculate base functions on flat prism.
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Partitioned (Indexed) Finite Element in Problem.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
SideNumber_multiIndex & getSideNumberTable() const
EntityHandle getEnt() const