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
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