v0.10.0
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 893 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpTractionBc()

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

Definition at line 894 of file EshelbianPlasticity.hpp.

895  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
896  bcFe(bc_fe) {}
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator

Member Function Documentation

◆ doWork()

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

Definition at line 893 of file EshelbianOperators.cpp.

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

Member Data Documentation

◆ bcFe

FeTractionBc& EshelbianPlasticity::OpTractionBc::bcFe
protected

Definition at line 900 of file EshelbianPlasticity.hpp.

◆ matPP

MatrixDouble EshelbianPlasticity::OpTractionBc::matPP
protected

Definition at line 901 of file EshelbianPlasticity.hpp.

◆ vecPv

MatrixDouble EshelbianPlasticity::OpTractionBc::vecPv
protected

Definition at line 902 of file EshelbianPlasticity.hpp.


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