529 {
531
535 }
536 if (row_data.getIndices().size() == 0)
538 if (col_data.getIndices().size() == 0)
540 int nb_dofs_row = row_data.getIndices().size();
541 int nb_dofs_col = col_data.getIndices().size();
544
545 kMatrix0.resize(nb_dofs_row, nb_dofs_col,
false);
546 kMatrix1.resize(nb_dofs_row, nb_dofs_col,
false);
547 kMatrix.resize(nb_dofs_row, nb_dofs_col,
false);
549
550 try {
551
552 int gg = 0;
553 for (int ff = 0; ff < 4; ff++) {
555 continue;
561 for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++, gg++) {
574 3, ublas::shallow_array_adaptor<double>(
576 double area = cblas_dnrm2(3, &normal[0], 1);
577
581
582
583 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
584 double n_row = row_data.getN()(gg, dd1);
585 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
586 double n_col = col_data.getN()(gg, dd2);
587 for (int dd3 = 0; dd3 < 3; dd3++) {
588 for (int dd4 = 0; dd4 < 3; dd4++) {
589 if (!P(dd3, dd4))
590 continue;
591 kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
592 val * n_row * P(dd3, dd4) * n_col * area;
593 }
594 }
595 }
596 }
597 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
598 double n_row = row_data.getN()(gg, dd1);
599 for (int dd2 = 0; dd2 < nb_dofs_col; dd2++) {
600 for (int dd3 = 0; dd3 < 3; dd3++) {
601 double t = cblas_ddot(3, &P(dd3, 0), 1, &
tRac_u(0, dd2),
603 kMatrix1(3 * dd1 + dd3, dd2) -= val * n_row *
t;
604 }
605 }
606 }
607 for (int dd1 = 0; dd1 < nb_dofs_row; dd1++) {
608 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
609 double n_col = col_data.getN()(gg, dd2);
610 for (int dd3 = 0; dd3 < 3; dd3++) {
611 double t = cblas_ddot(3, &P(0, dd3), 3, &
tRac_v(0, dd1),
614 }
615 }
616 }
617
621 (unsigned int)nb_face_gauss_pts) {
623 if (dP.size1() == 3 &&
624 dP.size2() == (unsigned int)nb_dofs_col) {
628 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
629 double n_row = row_data.getN()(gg, dd1);
630 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
631 for (int dd3 = 0; dd3 < 3; dd3++) {
632 for (int dd4 = 0; dd4 < 3; dd4++) {
633 kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
634 val * n_row * u[dd4] * dP(dd4, 3 * dd2 + dd4);
635 }
636 }
637 }
638 }
641 normal);
642 for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
643 double n_row = row_data.getN()(gg, dd1);
644 for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
645 for (int dd3 = 0; dd3 < 3; dd3++) {
646 for (int dd4 = 0; dd4 < 3; dd4++) {
647 kMatrix1(3 * dd1 + dd3, 3 * dd2 + dd4) +=
648 val * n_row *
t[dd4] * dP(dd4, 3 * dd2 + dd4);
649 }
650 }
651 }
652 }
653 }
654 }
655 }
656 }
657 }
658
661 }
662
665 "wrong number of gauss pts");
666 }
667
669 getFEMethod()->snes_B, nb_dofs_row, &row_data.getIndices()[0],
670 nb_dofs_col, &col_data.getIndices()[0], &
kMatrix(0, 0), ADD_VALUES);
672
673 } catch (const std::exception &ex) {
674 std::ostringstream ss;
675 ss << "throw in method: " << ex.what() << std::endl;
676 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
677 }
678
680 }
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
VectorShallowArrayAdaptor< double > VectorAdaptor
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
constexpr double t
plate stiffness
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
double gamma
Penalty term, see .
double phi
Nitsche method parameter, see .
std::vector< boost::shared_ptr< const NumeredEntFiniteElement > > facesFePtr
std::vector< MatrixDouble > faceNormals
std::map< EntityType, std::vector< std::vector< MatrixDouble > > > dP
derivative of projection matrix in respect DoFs This is EntityType, face, gauss point at face....
std::vector< std::vector< MatrixDouble > > P
projection matrix
std::vector< MatrixDouble > faceGaussPts
virtual MoFEMErrorCode getFaceRadius(int ff)
BlockData & nitscheBlockData
CommonData & nitscheCommonData
virtual MoFEMErrorCode getGammaH(double gamma, int gg)
NonlinearElasticElement::CommonData & commonData
virtual MoFEMErrorCode calculateP(int gg, int fgg, int ff)
NonlinearElasticElement::BlockData & dAta
MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff, MatrixDouble &jac, MatrixDouble &trac)
Range tEts
constrains elements in block set
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< MatrixDouble3by3 > sTress