v0.9.1
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 835 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpTractionBc()

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

Definition at line 836 of file EshelbianPlasticity.hpp.

837  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
838  bcFe(bc_fe) {}
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator

Member Function Documentation

◆ doWork()

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

Definition at line 882 of file EshelbianOperators.cpp.

882  {
884 
887 
888  auto t_normal = getFTensor1Normal();
889  const double nrm2 = sqrt(t_normal(i) * t_normal(i));
890  FTensor::Tensor1<double, 3> t_unit_normal;
891  t_unit_normal(i) = t_normal(i) / nrm2;
892  int nb_dofs = data.getFieldData().size();
893  int nb_integration_pts = data.getN().size1();
894  int nb_base_functions = data.getN().size2() / 3;
895  double ts_t = getFEMethod()->ts_t;
896 
897  auto integrate_matrix = [&]() {
899 
900  auto t_w = getFTensor0IntegrationWeight();
901  matPP.resize(nb_dofs / 3, nb_dofs / 3, false);
902  matPP.clear();
903 
904  auto t_row_base_fun = data.getFTensor1N<3>();
905  for (int gg = 0; gg != nb_integration_pts; ++gg) {
906 
907  int rr = 0;
908  for (; rr != nb_dofs / 3; ++rr) {
909  const double a = t_row_base_fun(i) * t_unit_normal(i);
910  auto t_col_base_fun = data.getFTensor1N<3>(gg, 0);
911  for (int cc = 0; cc != nb_dofs / 3; ++cc) {
912  const double b = t_col_base_fun(i) * t_unit_normal(i);
913  matPP(rr, cc) += t_w * a * b;
914  ++t_col_base_fun;
915  }
916  ++t_row_base_fun;
917  }
918 
919  for (; rr != nb_base_functions; ++rr)
920  ++t_row_base_fun;
921 
922  ++t_w;
923  }
924 
926  };
927 
928  auto integrate_rhs = [&](auto &bc) {
930 
931  auto t_w = getFTensor0IntegrationWeight();
932  vecPv.resize(3, nb_dofs / 3, false);
933  vecPv.clear();
934 
935  auto t_row_base_fun = data.getFTensor1N<3>();
936  double ts_t = getFEMethod()->ts_t;
937 
938  for (int gg = 0; gg != nb_integration_pts; ++gg) {
939  int rr = 0;
940  for (; rr != nb_dofs / 3; ++rr) {
941  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
942  for (int dd = 0; dd != 3; ++dd)
943  if (bc.flags[dd])
944  vecPv(dd, rr) += t * bc.vals[dd];
945  ++t_row_base_fun;
946  }
947  for (; rr != nb_base_functions; ++rr)
948  ++t_row_base_fun;
949  ++t_w;
950  }
952  };
953 
954  auto integrate_rhs_cook = [&](auto &bc) {
956 
957  vecPv.resize(3, nb_dofs / 3, false);
958  vecPv.clear();
959 
960  auto t_w = getFTensor0IntegrationWeight();
961  auto t_row_base_fun = data.getFTensor1N<3>();
962  auto t_coords = getFTensor1CoordsAtGaussPts();
963 
964  for (int gg = 0; gg != nb_integration_pts; ++gg) {
965 
966  auto calc_tau = [](double y) {
967  y -= 44;
968  y /= (60 - 44);
969  return -y * (y - 1) / 0.25;
970  };
971 
972  const double tau = calc_tau(t_coords(1));
973 
974  int rr = 0;
975  for (; rr != nb_dofs / 3; ++rr) {
976  const double t = ts_t * t_w * t_row_base_fun(i) * t_unit_normal(i);
977  for (int dd = 0; dd != 3; ++dd)
978  if (bc.flags[dd])
979  vecPv(dd, rr) += t * tau * bc.vals[dd];
980  ++t_row_base_fun;
981  }
982 
983  for (; rr != nb_base_functions; ++rr)
984  ++t_row_base_fun;
985  ++t_w;
986  ++t_coords;
987  }
989  };
990 
991  // get entity of face
992  EntityHandle fe_ent = getFEEntityHandle();
993  for (auto &bc : *(bcFe.bcData)) {
994  if (bc.faces.find(fe_ent) != bc.faces.end()) {
995 
996  int nb_dofs = data.getFieldData().size();
997  if (nb_dofs) {
998 
999  CHKERR integrate_matrix();
1000  if (bc.blockName.compare(20, 4, "COOK") == 0)
1001  CHKERR integrate_rhs_cook(bc);
1002  else
1003  CHKERR integrate_rhs(bc);
1004 
1005  const auto info =
1006  lapack_dposv('L', nb_dofs / 3, 3, &*matPP.data().begin(),
1007  nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1008  if (info > 0)
1009  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1010  "The leading minor of order %i is "
1011  "not positive; definite;\nthe "
1012  "solution could not be computed",
1013  info);
1014 
1015  for (int dd = 0; dd != 3; ++dd)
1016  if (bc.flags[dd])
1017  for (int rr = 0; rr != nb_dofs / 3; ++rr)
1018  data.getFieldDofs()[3 * rr + dd]->getFieldData() = vecPv(dd, rr);
1019 
1020  // int nb_integration_pts = data.getN().size1();
1021  // auto t_row_base_fun = data.getFTensor1N<3>();
1022  // auto t_coords = getFTensor1CoordsAtGaussPts();
1023  // for (int gg = 0; gg != nb_integration_pts; ++gg) {
1024 
1025  // FTensor::Tensor1<double, 3> t_v;
1026  // t_v(i) = 0;
1027 
1028  // int rr = 0;
1029  // for (; rr != nb_dofs / 3; ++rr) {
1030  // for (int dd = 0; dd != 3; ++dd) {
1031  // const double v = data.getFieldDofs()[3 * rr +
1032  // dd]->getFieldData(); t_v(dd) += t_row_base_fun(i) *
1033  // t_unit_normal(i) * v;
1034  // }
1035  // ++t_row_base_fun;
1036  // }
1037  // for (; rr != nb_base_functions; ++rr)
1038  // ++t_row_base_fun;
1039  // ++t_coords;
1040 
1041  // cerr << t_v << endl;
1042  // }
1043  }
1044  }
1045  }
1046 
1048 }
boost::shared_ptr< TractionBcVec > bcData
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:601
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26

Member Data Documentation

◆ bcFe

FeTractionBc& EshelbianPlasticity::OpTractionBc::bcFe
protected

Definition at line 842 of file EshelbianPlasticity.hpp.

◆ matPP

MatrixDouble EshelbianPlasticity::OpTractionBc::matPP
protected

Definition at line 843 of file EshelbianPlasticity.hpp.

◆ vecPv

MatrixDouble EshelbianPlasticity::OpTractionBc::vecPv
protected

Definition at line 844 of file EshelbianPlasticity.hpp.


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