v0.8.23
Public Member Functions | Protected Attributes | List of all members
EshelbianPlasticity::OpTractionBc Struct Reference

#include <users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp>

Inherits FaceUserDataOperator.

Collaboration diagram for EshelbianPlasticity::OpTractionBc:
[legend]

Public Member Functions

 OpTractionBc (std::string row_field, FeTractionBc &bc_fe)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 

Protected Attributes

FeTractionBcbcFe
 
MatrixDouble matPP
 
MatrixDouble vecPv
 

Detailed Description

Definition at line 907 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpTractionBc()

EshelbianPlasticity::OpTractionBc::OpTractionBc ( std::string  row_field,
FeTractionBc bc_fe 
)

Definition at line 908 of file EshelbianPlasticity.hpp.

909  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
910  bcFe(bc_fe) {}
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator

Member Function Documentation

◆ doWork()

MoFEMErrorCode EshelbianPlasticity::OpTractionBc::doWork ( int  side,
EntityType  type,
EntData data 
)

Definition at line 599 of file EshelbianOperators.cpp.

599  {
601 
604 
605  auto t_normal = getFTensor1Normal();
606  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
607  FTensor::Tensor1<double, 3> t_unit_normal;
608  t_unit_normal(i) = t_normal(i) / nrm2;
609  int nb_dofs = data.getFieldData().size();
610  int nb_integration_pts = data.getN().size1();
611  int nb_base_functions = data.getN().size2() / 3;
612  double ts_t = getFEMethod()->ts_t;
613 
614  auto integrate_matrix = [&]() {
616 
617  auto t_w = getFTensor0IntegrationWeight();
618  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
619  matPP.clear();
620 
621  auto t_row_base_fun = data.getFTensor1N<3>();
622  for (int gg = 0; gg != nb_integration_pts; ++gg) {
623 
624  int rr = 0;
625  for (; rr != nb_dofs / 3; ++rr) {
626  const double a = t_row_base_fun(i) * t_unit_normal(i);
627  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
628  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
629  const double b = t_col_base_fun(i) * t_unit_normal(i);
630  matPP(rr, cc) += t_w * a * b;
631  ++t_col_base_fun;
632  }
633  ++t_row_base_fun;
634  }
635 
636  for (; rr != nb_base_functions; ++rr)
637  ++t_row_base_fun;
638 
639  ++t_w;
640  }
641 
643  };
644 
645  auto integrate_rhs = [&](auto &bc) {
647 
648  auto t_w = getFTensor0IntegrationWeight();
649  vecPv.resize(3, nb_dofs / 3, false);
650  vecPv.clear();
651 
652  auto t_row_base_fun = data.getFTensor1N<3>();
653  double ts_t = getFEMethod()->ts_t;
654 
655  for (int gg = 0; gg != nb_integration_pts; ++gg) {
656  int rr = 0;
657  for (; rr != nb_dofs / 3; ++rr) {
658  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
659  for (int dd = 0; dd != 3; ++dd)
660  if (bc.flags[dd])
661  vecPv(dd, rr) += t * bc.vals[dd];
662  ++t_row_base_fun;
663  }
664  for (; rr != nb_base_functions; ++rr)
665  ++t_row_base_fun;
666  ++t_w;
667  }
669  };
670 
671  auto integrate_rhs_cook = [&](auto &bc) {
673 
674  vecPv.resize(3, nb_dofs / 3, false);
675  vecPv.clear();
676 
677  auto t_w = getFTensor0IntegrationWeight();
678  auto t_row_base_fun = data.getFTensor1N<3>();
679  auto t_coords = getFTensor1CoordsAtGaussPts();
680 
681  for (int gg = 0; gg != nb_integration_pts; ++gg) {
682 
683  auto calc_tau = [](double y) {
684  y -= 44;
685  y /= (60 - 44);
686  return -y * (y - 1) / 0.25;
687  };
688 
689  const double tau = calc_tau(t_coords(1));
690 
691  int rr = 0;
692  for (; rr != nb_dofs / 3; ++rr) {
693  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
694  for (int dd = 0; dd != 3; ++dd)
695  if (bc.flags[dd])
696  vecPv(dd, rr) += t * tau * bc.vals[dd];
697  ++t_row_base_fun;
698  }
699 
700  for (; rr != nb_base_functions; ++rr)
701  ++t_row_base_fun;
702  ++t_w;
703  ++t_coords;
704  }
706  };
707 
708  // get entity of face
709  EntityHandle fe_ent = getFEEntityHandle();
710  for (auto &bc : *(bcFe.bcData)) {
711  if (bc.faces.find(fe_ent) != bc.faces.end()) {
712 
713  int nb_dofs = data.getFieldData().size();
714  if (nb_dofs) {
715 
716  CHKERR integrate_matrix();
717  if (bc.blockName.compare(20, 4, "COOK") == 0)
718  CHKERR integrate_rhs_cook(bc);
719  else
720  CHKERR integrate_rhs(bc);
721 
722  const auto info =
723  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
724  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
725  if (info > 0)
726  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
727  "The leading minor of order %i is "
728  "not positive; definite;\nthe "
729  "solution could not be computed",
730  info);
731 
732  for (int dd = 0; dd != 3; ++dd)
733  if (bc.flags[dd])
734  for (int rr = 0; rr != nb_dofs / 3; ++rr)
735  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
736 
737  // int nb_integration_pts = data.getN().size1();
738  // auto t_row_base_fun = data.getFTensor1N<3>();
739  // auto t_coords = getFTensor1CoordsAtGaussPts();
740  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
741 
742  // FTensor::Tensor1<double, 3> t_v;
743  // t_v(i) = 0;
744 
745  // int rr = 0;
746  // for (; rr != nb_dofs / 3; ++rr) {
747  // for (int dd = 0; dd != 3; ++dd) {
748  // const double v = data.getFieldDofs()[3 * rr +
749  // dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
750  // t_unit_normal(i) * v;
751  // }
752  // ++t_row_base_fun;
753  // }
754  // for (; rr != nb_base_functions; ++rr)
755  // ++t_row_base_fun;
756  // ++t_coords;
757 
758  // cerr << t_v << endl;
759  // }
760  }
761  }
762  }
763 
765 }
boost::shared_ptr< TractionBcVec > bcData
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:223
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

Member Data Documentation

◆ bcFe

FeTractionBc& EshelbianPlasticity::OpTractionBc::bcFe
protected

Definition at line 914 of file EshelbianPlasticity.hpp.

◆ matPP

MatrixDouble EshelbianPlasticity::OpTractionBc::matPP
protected

Definition at line 915 of file EshelbianPlasticity.hpp.

◆ vecPv

MatrixDouble EshelbianPlasticity::OpTractionBc::vecPv
protected

Definition at line 916 of file EshelbianPlasticity.hpp.


The documentation for this struct was generated from the following files: