v0.13.0
ContactOps.hpp
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
/**
* \file ContactOps.hpp
* \example ContactOps.hpp
*/
namespace ContactOps {
//! [Common data]
struct CommonData {
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
boost::shared_ptr<MatrixDouble> contactStressPtr;
boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
boost::shared_ptr<MatrixDouble> contactTractionPtr;
boost::shared_ptr<MatrixDouble> contactDispPtr;
boost::shared_ptr<MatrixDouble> curlContactStressPtr;
};
//! [Common data]
struct OpConstrainBoundaryRhs : public AssemblyBoundaryEleOp {
OpConstrainBoundaryRhs(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpConstrainBoundaryLhs_dU : public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpConstrainBoundaryLhs_dTraction : public AssemblyBoundaryEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <typename T1, typename T2>
t_normal(i) = 0;
t_normal(1) = 1.;
return t_normal;
}
template <typename T>
inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
return (-0.5 - t_coords(1)) * t_normal(1);
}
template <typename T>
inline double gap(FTensor::Tensor1<T, SPACE_DIM> &t_disp,
return t_disp(i) * t_normal(i);
}
template <typename T>
return -t_traction(i) * t_normal(i);
}
inline double sign(double x) {
if (x == 0)
return 0;
else if (x > 0)
return 1;
else
return -1;
};
inline double w(const double g, const double t) { return g - cn * t; }
inline double constrian(double &&g0, double &&g, double &&t) {
return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
};
inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
return -cn * (1 + sign(w(g - g0, t))) / 2;
}
inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
return (1 + sign(w(g - g0, t))) / 2;
}
const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyBoundaryEleOp(field_name, field_name, DomainEleOp::OPROW),
commonDataPtr(common_data_ptr) {}
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_w = getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
size_t nb_base_functions = data.getN().size2() / 3;
auto t_base = data.getFTensor1N<3>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
t_rhs_constrains(i) =
t_contact_normal(i) *
constrian(gap0(t_coords, t_contact_normal),
gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
t_rhs_tangent_traction;
t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
t_rhs_tangent_traction(i) = cn * t_Q(i, j) * t_traction(j);
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++bb) {
const double beta = alpha * (t_base(i) * t_normal(i));
t_nf(i) -= beta * t_rhs_constrains(i);
t_nf(i) -= beta * t_rhs_tangent_disp(i);
t_nf(i) += beta * t_rhs_tangent_traction(i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
commonDataPtr(common_data_ptr) {
sYmm = false;
}
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
auto diff_constrain = diff_constrains_dgap(
gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_DIM * rr);
const double row_base = t_row_base(i) * t_normal(i);
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
++cc) {
const double beta = alpha * row_base * t_col_base;
t_mat(i, j) -= (beta * diff_constrain) * t_P(i, j);
t_mat(i, j) -= beta * t_Q(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
commonDataPtr(common_data_ptr) {
sYmm = false;
}
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
const double diff_traction = diff_constrains_dtraction(
gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_DIM * rr);
const double row_base = t_row_base(i) * t_normal(i);
auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
++cc) {
const double col_base = t_col_base(i) * t_normal(i);
const double beta = alpha * row_base * col_base;
t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
t_mat(i, j) += beta * cn * t_Q(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
}; // namespace ContactOps
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::OpBase AssemblyBoundaryEleOp
Definition: contact.cpp:73
constexpr double alpha
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:98
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:41
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:86
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:40
double w(const double g, const double t)
Definition: ContactOps.hpp:112
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:122
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:39
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:118
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:42
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:114
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:92
double sign(double x)
Definition: ContactOps.hpp:103
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
double cn
Definition: plastic.cpp:111
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr double g
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:35
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:25
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:26
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:27
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:32
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:33
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:28
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:30
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:31
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:277
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:72
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:286
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:61
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:199
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:208
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:126
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:50
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:132