40 int nb_gauss_pts = pts.size2();
42 auto ©_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
43 auto ©_diff_base_fun =
44 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
45 if (copy_base_fun.size1() != nb_gauss_pts)
47 "Inconsistent number of integration pts");
49 auto add_base_on_verts = [&] {
51 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts,
52 copy_base_fun.size2());
53 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(
54 nb_gauss_pts, copy_diff_base_fun.size2());
55 noalias(data.dataOnEntities[MBVERTEX][0].getN(base)) = copy_base_fun;
56 noalias(data.dataOnEntities[MBVERTEX][0].getDiffN(base)) =
62 auto add_base_on_edges = [&] {
64 std::array<int, 12> sense;
65 std::array<int, 12>
order;
66 if (data.spacesOnEntities[MBEDGE].test(
H1)) {
67 if (data.dataOnEntities[MBEDGE].size() != 12)
69 "Expected 12 data on entities");
71 std::array<double *, 12> h1_edge_n;
72 std::array<double *, 12> diff_h1_egde_n;
73 bool nb_dofs_sum =
false;
74 for (
int ee = 0; ee != 12; ++ee) {
75 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
77 "Sense of entity not set");
79 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
80 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
82 int nb_dofs =
NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
83 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
85 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
86 nb_gauss_pts, 3 * nb_dofs,
false);
88 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
90 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
92 nb_dofs_sum |= (nb_dofs > 0);
96 sense.data(),
order.data(), &*copy_base_fun.data().begin(),
97 &*copy_diff_base_fun.data().begin(), h1_edge_n.data(),
98 diff_h1_egde_n.data(), nb_gauss_pts);
101 for (
int ee = 0; ee != 12; ++ee) {
102 data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0,
false);
103 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0,
false);
110 auto add_base_on_quads = [&]() {
112 std::array<int, 6>
order;
113 if (data.spacesOnEntities[MBQUAD].test(
H1)) {
115 if (data.dataOnEntities[MBQUAD].size() != 6)
117 "Expected six faces on hex");
119 std::array<double *, 6> h1_face_n;
120 std::array<double *, 6> diff_h1_face_n;
121 bool nb_dofs_sum =
false;
122 for (
int ff = 0; ff != 6; ++ff) {
124 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
127 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
129 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(
130 nb_gauss_pts, 3 * nb_dofs,
false);
133 &*data.dataOnEntities[MBQUAD][ff].getN(base).data().begin();
135 &*data.dataOnEntities[MBQUAD][ff].getDiffN(base).data().begin();
137 nb_dofs_sum |= (nb_dofs > 0);
139 if (data.facesNodes.size1() != 6)
141 "Expected six face nodes");
142 if (data.facesNodes.size2() != 4)
144 "Expected four nodes on face");
148 &*data.facesNodes.data().begin(),
149 &*data.facesNodesOrder.data().begin(),
order.data(),
150 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
151 h1_face_n.data(), diff_h1_face_n.data(), nb_gauss_pts);
155 for (
int ff = 0; ff != 6; ++ff) {
156 data.dataOnEntities[MBQUAD][ff].getN(base).resize(0,
false);
157 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(0, 0,
false);
165 auto add_base_on_volume = [&]() {
168 if (data.spacesOnEntities[MBHEX].test(
H1)) {
170 int order = data.dataOnEntities[MBHEX][0].getOrder();
172 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
174 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(
175 nb_gauss_pts, 3 * nb_vol_dofs,
false);
180 p.data(), &*copy_base_fun.data().begin(),
181 &*copy_diff_base_fun.data().begin(),
182 &*data.dataOnEntities[MBHEX][0].getN(base).data().begin(),
183 &*data.dataOnEntities[MBHEX][0].getDiffN(base).data().begin(),
187 data.dataOnEntities[MBHEX][0].getN(base).resize(0, 0,
false);
188 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(0, 0,
false);
194 CHKERR add_base_on_verts();
195 CHKERR add_base_on_edges();
196 CHKERR add_base_on_quads();
197 CHKERR add_base_on_volume();
223 int nb_gauss_pts = pts.size2();
225 auto ©_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
226 auto ©_diff_base_fun =
227 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
229 if (nb_gauss_pts != copy_base_fun.size1())
231 "Wrong number of gauss pts");
233 if (nb_gauss_pts != copy_diff_base_fun.size1())
235 "Wrong number of gauss pts");
237 auto &base_fun = data.dataOnEntities[MBHEX][0].getN(base);
238 auto &diff_base_fun = data.dataOnEntities[MBHEX][0].getDiffN(base);
239 const int vol_order = data.dataOnEntities[MBHEX][0].getOrder();
242 base_fun.resize(nb_gauss_pts, nb_dofs,
false);
243 diff_base_fun.resize(nb_gauss_pts, 3 * nb_dofs,
false);
246 const std::array<int, 3>
p{vol_order, vol_order, vol_order};
248 p.data(), &*copy_base_fun.data().begin(),
249 &*copy_diff_base_fun.data().begin(), &*base_fun.data().begin(),
250 &*diff_base_fun.data().begin(), nb_gauss_pts);
276 const int nb_gauss_pts = pts.size2();
278 auto ©_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
279 auto ©_diff_base_fun =
280 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
282 if (nb_gauss_pts != copy_base_fun.size1())
284 "Wrong number of gauss pts");
286 if (nb_gauss_pts != copy_diff_base_fun.size1())
288 "Wrong number of gauss pts");
291 if (data.spacesOnEntities[MBQUAD].test(
HDIV)) {
293 if (data.dataOnEntities[MBQUAD].size() != 6)
295 "Expected six data structures on Hex");
297 std::array<int, 6>
order;
298 std::array<double *, 6> hdiv_face_n;
299 std::array<double *, 6> diff_hdiv_face_n;
301 bool sum_nb_dofs =
false;
302 for (
int ff = 0; ff != 6; ff++) {
303 if (data.dataOnEntities[MBQUAD][ff].getSense() == 0)
305 "Sense pn quad <%d> not set", ff);
307 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
308 if (data.facesNodes.size1() != 6)
310 "Expected six faces");
311 if (data.facesNodes.size2() != 4)
313 "Expected four nodes on face");
316 auto &face_n = data.dataOnEntities[MBQUAD][ff].getN(base);
317 auto &diff_face_n = data.dataOnEntities[MBQUAD][ff].getDiffN(base);
318 face_n.resize(nb_gauss_pts, 3 * nb_dofs,
false);
319 diff_face_n.resize(nb_gauss_pts, 9 * nb_dofs,
false);
321 hdiv_face_n[ff] = &*face_n.data().begin();
322 diff_hdiv_face_n[ff] = &*diff_face_n.data().begin();
324 sum_nb_dofs |= (nb_dofs > 0);
329 &*data.facesNodes.data().begin(),
330 &*data.facesNodesOrder.data().begin(),
order.data(),
331 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
332 hdiv_face_n.data(), diff_hdiv_face_n.data(), nb_gauss_pts);
335 for (
int ff = 0; ff != 6; ff++) {
336 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
338 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
344 for (
int ff = 0; ff != 6; ff++) {
345 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
false);
346 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
352 if (data.spacesOnEntities[MBHEX].test(
HDIV)) {
354 const int order = data.dataOnEntities[MBHEX][0].getOrder();
355 const int nb_dofs_family =
358 volFamily.resize(3, 3 * nb_dofs_family * nb_gauss_pts);
360 if (nb_dofs_family) {
364 std::array<double *, 3> diff_family_ptr = {
369 p.data(), &*copy_base_fun.data().begin(),
370 &*copy_diff_base_fun.data().begin(), family_ptr.data(),
371 diff_family_ptr.data(), nb_gauss_pts);
374 auto &face_n = data.dataOnEntities[MBHEX][0].getN(base);
375 auto &diff_face_n = data.dataOnEntities[MBHEX][0].getDiffN(base);
376 face_n.resize(nb_gauss_pts, 3 * nb_vol_dofs,
false);
377 diff_face_n.resize(nb_gauss_pts, 9 * nb_vol_dofs,
false);
382 double *ptr = &face_n(0, 0);
384 for (
int j = 0;
j != 3; ++
j) {
389 for (
int j = 0;
j != 3; ++
j) {
394 for (
int j = 0;
j != 3; ++
j) {
404 double *diff_ptr = &diff_face_n(0, 0);
406 for (
int j = 0;
j != 9; ++
j) {
407 *diff_ptr = *diff_ptr_f0;
411 for (
int j = 0;
j != 9; ++
j) {
412 *diff_ptr = *diff_ptr_f1;
416 for (
int j = 0;
j != 9; ++
j) {
417 *diff_ptr = *diff_ptr_f2;
423 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0,
false);
424 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
429 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0,
false);
430 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
false);
457 const int nb_gauss_pts = pts.size2();
459 auto ©_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
460 auto ©_diff_base_fun =
461 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
463 if (nb_gauss_pts != copy_base_fun.size1())
465 "Wrong number of gauss pts");
467 if (nb_gauss_pts != copy_diff_base_fun.size1())
469 "Wrong number of gauss pts");
472 if (data.spacesOnEntities[MBEDGE].test(
HCURL)) {
473 std::array<int, 12> sense;
474 std::array<int, 12>
order;
475 if (data.dataOnEntities[MBEDGE].size() != 12)
477 "Expected 12 edges data structures on Hex");
479 std::array<double *, 12> hcurl_edge_n;
480 std::array<double *, 12> diff_hcurl_edge_n;
481 bool sum_nb_dofs =
false;
483 for (
int ee = 0; ee != 12; ++ee) {
484 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
486 "Sense on edge <%d> on Hex not set", ee);
488 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
489 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
491 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
493 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
496 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
497 diff_hcurl_edge_n[ee] =
498 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
500 sum_nb_dofs |= (nb_dofs > 0);
505 sense.data(),
order.data(), &*copy_base_fun.data().begin(),
506 &*copy_diff_base_fun.data().begin(), hcurl_edge_n.data(),
507 diff_hcurl_edge_n.data(), nb_gauss_pts);
511 for (
int ee = 0; ee != 12; ee++) {
512 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0,
false);
513 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
519 if (data.spacesOnEntities[MBQUAD].test(
HCURL)) {
521 if (data.dataOnEntities[MBQUAD].size() != 6)
523 "Expected six data structures on Hex");
525 std::array<int, 6>
order;
526 double *face_family_ptr[6][2];
527 double *diff_face_family_ptr[6][2];
529 bool sum_nb_dofs =
false;
530 for (
int ff = 0; ff != 6; ff++) {
531 if (data.dataOnEntities[MBQUAD][ff].getSense() == 0)
533 "Sense pn quad <%d> not set", ff);
535 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
536 if (data.facesNodes.size1() != 6)
538 "Expected six faces");
539 if (data.facesNodes.size2() != 4)
541 "Expected four nodes on face");
543 const int nb_family_dofs =
545 faceFamily[ff].resize(2, 3 * nb_family_dofs * nb_gauss_pts,
false);
546 diffFaceFamily[ff].resize(2, 9 * nb_family_dofs * nb_gauss_pts,
false);
548 if (nb_family_dofs) {
549 face_family_ptr[ff][0] = &(
faceFamily[ff](0, 0));
550 face_family_ptr[ff][1] = &(
faceFamily[ff](1, 0));
555 sum_nb_dofs |= (nb_family_dofs > 0);
560 &*data.facesNodes.data().begin(),
561 &*data.facesNodesOrder.data().begin(),
order.data(),
562 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
563 face_family_ptr, diff_face_family_ptr, nb_gauss_pts);
565 for (
int ff = 0; ff != 6; ++ff) {
571 auto &face_n = data.dataOnEntities[MBQUAD][ff].getN(base);
572 auto &diff_face_n = data.dataOnEntities[MBQUAD][ff].getDiffN(base);
573 face_n.resize(nb_gauss_pts, 3 * nb_dofs,
false);
574 diff_face_n.resize(nb_gauss_pts, 9 * nb_dofs,
false);
578 double *ptr = &face_n(0, 0);
580 for (
int j = 0;
j != 3; ++
j) {
585 for (
int j = 0;
j != 3; ++
j) {
594 double *diff_ptr = &diff_face_n(0, 0);
596 for (
int j = 0;
j != 9; ++
j) {
597 *diff_ptr = *diff_ptr_f0;
601 for (
int j = 0;
j != 9; ++
j) {
602 *diff_ptr = *diff_ptr_f1;
610 for (
int ff = 0; ff != 6; ff++) {
611 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
613 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
619 for (
int ff = 0; ff != 6; ff++) {
620 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
false);
621 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
627 if (data.spacesOnEntities[MBHEX].test(
HCURL)) {
629 const int order = data.dataOnEntities[MBHEX][0].getOrder();
632 volFamily.resize(3, 3 * nb_dofs * nb_gauss_pts);
638 std::array<double *, 3> diff_family_ptr = {
643 p.data(), &*copy_base_fun.data().begin(),
644 &*copy_diff_base_fun.data().begin(), family_ptr.data(),
645 diff_family_ptr.data(), nb_gauss_pts);
648 auto &face_n = data.dataOnEntities[MBHEX][0].getN(base);
649 auto &diff_face_n = data.dataOnEntities[MBHEX][0].getDiffN(base);
650 face_n.resize(nb_gauss_pts, 3 * nb_vol_dofs,
false);
651 diff_face_n.resize(nb_gauss_pts, 9 * nb_vol_dofs,
false);
656 double *ptr = &face_n(0, 0);
658 for (
int j = 0;
j != 3; ++
j) {
663 for (
int j = 0;
j != 3; ++
j) {
668 for (
int j = 0;
j != 3; ++
j) {
678 double *diff_ptr = &diff_face_n(0, 0);
680 for (
int j = 0;
j != 9; ++
j) {
681 *diff_ptr = *diff_ptr_f0;
685 for (
int j = 0;
j != 9; ++
j) {
686 *diff_ptr = *diff_ptr_f1;
690 for (
int j = 0;
j != 9; ++
j) {
691 *diff_ptr = *diff_ptr_f2;
697 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0,
false);
698 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
703 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0,
false);
704 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
false);
712 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
717 int nb_gauss_pts = pts.size2();
724 "Wrong dimension of pts, should be at least 3 rows with coordinates");
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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
#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 ...
FTensor::Index< 'n', SPACE_DIM > n
#define NBFACEQUAD_DEMKOWICZ_HDIV(P)
#define NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL(P, Q)
Number of base functions on quad for Hcurl space.
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBVOLUMEHEX_DEMKOWICZ_HCURL(P)
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
#define NBVOLUMEHEX_H1(P)
Number of base functions on hex for H1 space.
#define NBVOLUMEHEX_DEMKOWICZ_HDIV(P)
#define NBFACEQUAD_DEMKOWICZ_HCURL(P)
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBVOLUMEHEX_L2(P)
Number of base functions on hexahedron for L2 space.
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, Q, R)
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, Q, R)
FTensor::Index< 'j', 3 > j
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
MoFEMErrorCode H1_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
MoFEMErrorCode L2_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *volN, double *diff_volN, int nb_integration_pts)
MoFEMErrorCode H1_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6], double *diff_faceN[6], int nb_integration_pts)
MoFEMErrorCode Hdiv_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN, double *faceN[6], double *div_faceN[6], int nb_integration_pts)
MoFEMErrorCode Hcurl_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6][2], double *diff_faceN[6][2], int nb_integration_pts)
MoFEMErrorCode H1_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN, int nb_integration_pts)
MoFEMErrorCode Hdiv_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *bubleN[3], double *div_bubleN[3], int nb_integration_pts)
MoFEMErrorCode Hcurl_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *volN[3], double *diff_volN[3], int nb_integration_pts)
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
const FieldApproximationBase bAse
Calculate base functions on tetrahedral.
MoFEMErrorCode getValueH1DemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MatrixDouble diffVolFamily
std::array< MatrixDouble, 6 > faceFamily
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
MoFEMErrorCode getValueL2DemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
std::array< MatrixDouble, 6 > diffFaceFamily
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
EntPolynomialBaseCtx * cTx
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.