v0.14.0
Public Member Functions | Public Attributes | List of all members
FractureMechanics::OpAnalyticalMaterialTraction Struct Reference

#include <users_modules/fracture_mechanics/src/CrackFrontElement.hpp>

Inheritance diagram for FractureMechanics::OpAnalyticalMaterialTraction:
[legend]
Collaboration diagram for FractureMechanics::OpAnalyticalMaterialTraction:
[legend]

Public Member Functions

 OpAnalyticalMaterialTraction (MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity, const std::string field_name, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &analytical_force_op, Range *forces_only_on_entities_row=NULL)
 
MoFEMErrorCode doWork (int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
 

Public Attributes

const Range tRis
 
boost::ptr_vector< MethodForForceScaling > & methodsOp
 
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
 
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
 
boost::shared_ptr< MatrixDouble > hG
 spatial deformation gradient More...
 
boost::shared_ptr< MatrixDouble > HG
 spatial deformation gradient More...
 
VectorDouble sNrm
 Length of the normal vector. More...
 
bool setSingularCoordinates
 
Range forcesOnlyOnEntitiesRow
 
VectorInt iNdices
 
VectorDouble Nf
 

Detailed Description

Definition at line 444 of file CrackFrontElement.hpp.

Constructor & Destructor Documentation

◆ OpAnalyticalMaterialTraction()

FractureMechanics::OpAnalyticalMaterialTraction::OpAnalyticalMaterialTraction ( MoFEM::Interface &  m_field,
const bool &  set_singular_coordinates,
const Range &  crack_front_nodes,
const Range &  crack_front_nodes_edges,
const Range &  crack_front_elements,
PetscBool  add_singularity,
const std::string  field_name,
const Range  tris,
boost::ptr_vector< MethodForForceScaling > &  methods_op,
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &  analytical_force_op,
Range *  forces_only_on_entities_row = NULL 
)

Definition at line 992 of file CrackFrontElement.cpp.

1001  : FaceElementForcesAndSourcesCore::UserDataOperator(
1002  field_name, UserDataOperator::OPROW),
1003  tRis(tris), methodsOp(methods_op), analyticalForceOp(analytical_force_op),
1004  volSideFe(m_field, set_singular_coordinates, crack_front_nodes,
1005  crack_front_nodes_edges, crack_front_elements, add_singularity),
1006  setSingularCoordinates(set_singular_coordinates) {
1007  hG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1008  HG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1009  volSideFe.getOpPtrVector().push_back(
1010  new OpCalculateVectorFieldGradient<3, 3>("SPATIAL_POSITION", hG));
1011  volSideFe.getOpPtrVector().push_back(
1012  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS", HG));
1013  volSideFe.meshPositionsFieldName = "NONE";
1014  if (forces_only_on_entities_row) {
1015  forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
1016  }
1017 }

Member Function Documentation

◆ doWork()

MoFEMErrorCode FractureMechanics::OpAnalyticalMaterialTraction::doWork ( int  side,
EntityType  type,
DataForcesAndSourcesCore::EntData &  data 
)

Definition at line 1020 of file CrackFrontElement.cpp.

1021  {
1022  MoFEMFunctionBeginHot;
1023 
1024  if (data.getIndices().size() == 0)
1025  MoFEMFunctionReturnHot(0);
1026  EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
1027  if (tRis.find(ent) == tRis.end())
1028  MoFEMFunctionReturnHot(0);
1029 
1030  if (type == MBVERTEX) {
1031  CHKERR loopSideVolumes("MATERIAL", volSideFe);
1032  }
1033  if (type != MBVERTEX) {
1034  if (volSideFe.singularElement) {
1035  const EntityHandle ent = data.getFieldDofs()[0]->getEnt();
1036  int side_number, sense, offset;
1037  CHKERR volSideFe.mField.get_moab().side_number(
1038  volSideFe.numeredEntFiniteElementPtr->getEnt(), ent, side_number,
1039  sense, offset);
1040  data.getN() =
1041  volSideFe.getEntData(H1, type, side_number).getN(data.getBase());
1042  }
1043  } else if (volSideFe.singularElement) {
1044  // Set value of base functions at singular integration points
1045  // FIXME: That need to be moved to separate operator such that values
1046  // of base function are appropriately evaluated
1047  for (unsigned int gg = 0; gg != volSideFe.gaussPts.size2(); gg++) {
1048  for (int dd = 0; dd != 3; dd++) {
1049  volSideFe.gaussPts(dd, gg) += volSideFe.singularRefCoords(gg, dd);
1050  }
1051  }
1052  CHKERR ShapeMBTET(
1053  &*volSideFe.getEntData(H1, MBVERTEX, 0).getN(NOBASE).data().begin(),
1054  &volSideFe.gaussPts(0, 0), &volSideFe.gaussPts(1, 0),
1055  &volSideFe.gaussPts(2, 0), volSideFe.gaussPts.size2());
1056  CHKERR volSideFe.calHierarchicalBaseFunctionsOnElement(data.getBase());
1057  map<int, int> map_face_nodes;
1058  for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
1059  for (int nn = 0; nn != 3; nn++) {
1060  data.getN()(gg, nn) =
1061  volSideFe.getEntData(H1, MBVERTEX, 0)
1062  .getN(data.getBase())(gg, volSideFe.getFaceConnMap()[nn]);
1063  }
1064  }
1065  }
1066 
1067  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1068  int nb_row_base = data.getIndices().size() / rank;
1069 
1070  Nf.resize(data.getIndices().size(), false);
1071  Nf.clear();
1072 
1073  VectorDouble3 coords(3);
1074  VectorDouble3 normal(3);
1075  VectorDouble3 force(3);
1076  VectorDouble3 FTforce(3);
1077 
1078  FTensor::Index<'i', 3> i;
1079  FTensor::Index<'j', 3> j;
1080  FTensor::Index<'k', 3> k;
1081 
1082  // If singular element apply transformation to h tensor
1083  if (type == MBVERTEX && volSideFe.singularElement) {
1084  // cerr << volSideFe.sJac << endl;
1085  // cerr << volSideFe.invSJac << endl;
1086  double *inv_jac_ptr = &*volSideFe.invSJac.data().begin();
1087  FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
1088  inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
1089  &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
1090  &inv_jac_ptr[8], 9);
1091  auto t_hg = getFTensor2FromMat<3, 3>(*hG);
1092  auto t_normal = getFTensor1Normal();
1093  double nrm = 2 * getArea();
1094  FTensor::Tensor1<double, 3> t_s_normal;
1095  sNrm.resize(data.getN().size1(), false);
1096  for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
1097  // Push normal to transformed singular element
1098  double det = volSideFe.detS[gg];
1099  t_s_normal(i) = det * t_inv_s_jac(j, i) * t_normal(j) / nrm;
1100  sNrm[gg] = sqrt(t_s_normal(i) * t_s_normal(i));
1101  FTensor::Tensor2<double, 3, 3> t_tmp;
1102  t_tmp(i, j) = t_hg(i, k) * t_inv_s_jac(k, j);
1103  t_hg(i, j) = t_tmp(i, j);
1104  ++t_inv_s_jac;
1105  ++t_hg;
1106  }
1107  }
1108 
1109  FTensor::Tensor1<double *, 3> t_force(&force[0], &force[1], &force[2]);
1110  FTensor::Tensor1<double *, 3> t_ft_force(&FTforce[0], &FTforce[1],
1111  &FTforce[2]);
1112  auto t_hg = getFTensor2FromMat<3, 3>(*hG);
1113  auto t_Hg = getFTensor2FromMat<3, 3>(*HG);
1114  FTensor::Tensor2<double, 3, 3> t_inv_Hg;
1115  FTensor::Tensor2<double, 3, 3> t_F;
1116 
1117  // cerr << "\n\n";
1118  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
1119 
1120  // get integration weight and Jacobian of integration point (area of face)
1121  double val = getGaussPts()(2, gg);
1122  val *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
1123  for (int dd = 0; dd != 3; dd++) {
1124  coords[dd] = getCoordsAtGaussPts()(gg, dd);
1125  normal[dd] = getNormalsAtGaussPts()(gg, dd);
1126  }
1127  if (volSideFe.singularElement) {
1128  for (int dd = 0; dd != 3; dd++) {
1129  coords[dd] += volSideFe.singularDisp(gg, dd);
1130  }
1131  val *= sNrm[gg];
1132  }
1133 
1134  for (boost::ptr_vector<
1135  NeumannForcesSurface::MethodForAnalyticalForce>::iterator vit =
1136  analyticalForceOp.begin();
1137  vit != analyticalForceOp.end(); vit++) {
1138 
1139  // Calculate force
1140  CHKERR vit->getForce(ent, coords, normal, force);
1141  // Invert matrix H
1142  double det;
1143  CHKERR determinantTensor3by3(t_Hg, det);
1144  CHKERR invertTensor3by3(t_Hg, det, t_inv_Hg);
1145 
1146  // Calculate gradient of deformation
1147  t_F(i, k) = t_hg(i, j) * t_inv_Hg(j, k);
1148  t_ft_force(i) = -t_F(j, i) * t_force(j);
1149 
1150  for (int rr = 0; rr != 3; rr++) {
1151  cblas_daxpy(nb_row_base, val * FTforce[rr], &data.getN()(gg, 0), 1,
1152  &Nf[rr], rank);
1153  }
1154  }
1155 
1156  ++t_hg;
1157  ++t_Hg;
1158  }
1159 
1160  // Scale force using user defined scaling operator
1161  CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
1162 
1163  // Assemble force into vector
1164  int nb_dofs = data.getIndices().size();
1165  int *indices_ptr = &data.getIndices()[0];
1166  if (!forcesOnlyOnEntitiesRow.empty()) {
1167  iNdices.resize(nb_dofs, false);
1168  noalias(iNdices) = data.getIndices();
1169  indices_ptr = &iNdices[0];
1170  auto &dofs = data.getFieldDofs();
1171  auto dit = dofs.begin();
1172  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1173  if (auto dof = (*dit)) {
1174  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
1175  forcesOnlyOnEntitiesRow.end()) {
1176  iNdices[ii] = -1;
1177  }
1178  }
1179  }
1180  }
1181  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &Nf[0],
1182  ADD_VALUES);
1183  MoFEMFunctionReturnHot(0);
1184 }

Member Data Documentation

◆ analyticalForceOp

boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>& FractureMechanics::OpAnalyticalMaterialTraction::analyticalForceOp

Definition at line 450 of file CrackFrontElement.hpp.

◆ forcesOnlyOnEntitiesRow

Range FractureMechanics::OpAnalyticalMaterialTraction::forcesOnlyOnEntitiesRow

Definition at line 458 of file CrackFrontElement.hpp.

◆ hG

boost::shared_ptr<MatrixDouble> FractureMechanics::OpAnalyticalMaterialTraction::hG

spatial deformation gradient

Definition at line 454 of file CrackFrontElement.hpp.

◆ HG

boost::shared_ptr<MatrixDouble> FractureMechanics::OpAnalyticalMaterialTraction::HG

spatial deformation gradient

Definition at line 455 of file CrackFrontElement.hpp.

◆ iNdices

VectorInt FractureMechanics::OpAnalyticalMaterialTraction::iNdices

Definition at line 460 of file CrackFrontElement.hpp.

◆ methodsOp

boost::ptr_vector<MethodForForceScaling>& FractureMechanics::OpAnalyticalMaterialTraction::methodsOp

Definition at line 448 of file CrackFrontElement.hpp.

◆ Nf

VectorDouble FractureMechanics::OpAnalyticalMaterialTraction::Nf

Definition at line 472 of file CrackFrontElement.hpp.

◆ setSingularCoordinates

bool FractureMechanics::OpAnalyticalMaterialTraction::setSingularCoordinates

Definition at line 457 of file CrackFrontElement.hpp.

◆ sNrm

VectorDouble FractureMechanics::OpAnalyticalMaterialTraction::sNrm

Length of the normal vector.

Definition at line 456 of file CrackFrontElement.hpp.

◆ tRis

const Range FractureMechanics::OpAnalyticalMaterialTraction::tRis

Definition at line 447 of file CrackFrontElement.hpp.

◆ volSideFe

CrackFrontSingularBase<FirendVolumeOnSide, VolumeElementForcesAndSourcesCore> FractureMechanics::OpAnalyticalMaterialTraction::volSideFe

Definition at line 453 of file CrackFrontElement.hpp.

The documentation for this struct was generated from the following files:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FractureMechanics::OpAnalyticalMaterialTraction::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: CrackFrontElement.hpp:458
H1
@ H1
continuous field
Definition: definitions.h:85
FractureMechanics::OpAnalyticalMaterialTraction::hG
boost::shared_ptr< MatrixDouble > hG
spatial deformation gradient
Definition: CrackFrontElement.hpp:454
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1< double, 3 >
EntityHandle
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::OpAnalyticalMaterialTraction::analyticalForceOp
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
Definition: CrackFrontElement.hpp:450
FractureMechanics::OpAnalyticalMaterialTraction::volSideFe
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
Definition: CrackFrontElement.hpp:453
FractureMechanics::OpAnalyticalMaterialTraction::Nf
VectorDouble Nf
Definition: CrackFrontElement.hpp:472
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
FTensor::Tensor2< double *, 3, 3 >
FractureMechanics::OpAnalyticalMaterialTraction::setSingularCoordinates
bool setSingularCoordinates
Definition: CrackFrontElement.hpp:457
FractureMechanics::OpAnalyticalMaterialTraction::HG
boost::shared_ptr< MatrixDouble > HG
spatial deformation gradient
Definition: CrackFrontElement.hpp:455
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1588
FractureMechanics::OpAnalyticalMaterialTraction::iNdices
VectorInt iNdices
Definition: CrackFrontElement.hpp:460
NeumannForcesSurface::MethodForAnalyticalForce
Analytical force method.
Definition: SurfacePressure.hpp:21
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FractureMechanics::OpAnalyticalMaterialTraction::sNrm
VectorDouble sNrm
Length of the normal vector.
Definition: CrackFrontElement.hpp:456
FractureMechanics::OpAnalyticalMaterialTraction::tRis
const Range tRis
Definition: CrackFrontElement.hpp:447
convert.type
type
Definition: convert.py:64
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
FTensor::dd
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
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FractureMechanics::OpAnalyticalMaterialTraction::methodsOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
Definition: CrackFrontElement.hpp:448
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
Generated by Doxygen 1.8.17 and hosted at University of Glasgow