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

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

Inherits FaceElementForcesAndSourcesCore::UserDataOperator.

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.

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}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
constexpr auto field_name
boost::shared_ptr< MatrixDouble > hG
spatial deformation gradient
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
boost::shared_ptr< MatrixDouble > HG
spatial deformation gradient
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.

Member Function Documentation

◆ doWork()

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

Definition at line 1020 of file CrackFrontElement.cpp.

1021 {
1023
1024 if (data.getIndices().size() == 0)
1026 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
1027 if (tRis.find(ent) == tRis.end())
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 }
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);
1075 VectorDouble3 force(3);
1076 VectorDouble3 FTforce(3);
1077
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();
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));
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);
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<
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
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()) ==
1176 iNdices[ii] = -1;
1177 }
1178 }
1179 }
1180 }
1181 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &Nf[0],
1182 ADD_VALUES);
1184}
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
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
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
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:1210
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
VectorDouble sNrm
Length of the normal vector.
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)

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: